#####READ BEFORE RUNNING CODE#####

#This code was developed using R version 4.2.1 for Mac
#It relies on the Zelig package, which is no longer updated or available on CRAN
#You must install Zelig 5.1.5 from the CRAN archive; this requires Matrix 1.2-14, which is also available on the CRAN archive
#Ensure that you remove existing versions of R and Matrix before installing the correct versions
#Libraries must be installed once before running the code, and can simply be loaded every time the file is reopened thereafter

#####INSTALL LIBRARIES#####

install.packages("Amelia")
install.packages("dplyr")
install.packages("MatchIt")
install.packages("MCMCpack")
install.packages("VGAM")
install.packages("survival") #Select no if prompted
install.packages("https://cran.r-project.org/src/contrib/Archive/Zelig/Zelig_5.1.5.tar.gz", repos = NULL, type = "source")
install.packages("gplots")
install.packages("Matching")
install.packages("weights")



#####LOAD LIBRARIES#####

rm(list=ls())
library(Zelig)
library(gplots)
library(Matching)
library(weights)


#####IMPORT DATA FILES#####

#Set working directory before proceeding

stacked <- read.csv("StackedAllConjoint.csv")
head(stacked)

stacked.w <- read.csv("StackedWhiteConjoint.csv") #Read in stacked dataset for whites
head(stacked.w)
nrow(stacked.w)

stacked.b <- read.csv("StackedBlackConjoint.csv") #Read in stacked dataset for Blacks
head(stacked.b)
nrow(stacked.b)


standard <- read.csv("StandardAll.csv") #Combine black and white standard datasets
head(standard)

standard.w <- read.csv("StandardWhite.csv") #Read in standard dataset for whites
head(standard.w)
nrow(standard.w)

standard.b <- read.csv("StandardBlack.csv") #Read in standard dataset for Blacks
head(standard.b)
nrow(standard.b)




#####MANIPULATION CHECK (FIGURE 2)#####

##Subset data##

cont.white <- subset(standard, treatment == 0 & race == "white") #Subset standard dataset by treatment condition and race
treat.white <- subset(standard, treatment == 1 & race == "white")
cont.black <- subset(standard, treatment == 0 & race == "black") 
treat.black <- subset(standard, treatment == 1 & race == "black")


cont.white.men <- subset(standard, treatment == 0 & female == 0 & race == "white") #Subset standard dataset by treatment, gender, and race
treat.white.men <- subset(standard, treatment == 1 & female == 0 & race == "white")
cont.white.women <- subset(standard, treatment == 0 & female == 1 & race == "white")
treat.white.women <- subset(standard, treatment == 1 & female == 1 & race == "white")
cont.black.men <- subset(standard, treatment == 0 & female == 0 & race == "black")
treat.black.men <- subset(standard, treatment == 1 & female == 0 & race == "black")
cont.black.women <- subset(standard, treatment == 0 & female == 1 & race == "black")
treat.black.women <- subset(standard, treatment == 1 & female == 1 & race == "black")


##Get summary statistics and t-tests##

#Compare whites in treatment and control
cont.white.mean <- mean(cont.white$clintonid, na.rm = TRUE) #Get means and standard errors for whites in control
cont.white.se <- sd(cont.white$clintonid, na.rm = TRUE)/sqrt(length(cont.white$clintonid))
print(cont.white.mean) #0.195
print(cont.white.se) #0.019

treat.white.mean <- mean(treat.white$clintonid, na.rm = TRUE) #Get means and standard errors for whites in treatment
treat.white.se <- sd(treat.white$clintonid, na.rm = TRUE)/sqrt(length(treat.white$clintonid))
print(treat.white.mean) #0.275
print(treat.white.se) #0.022

white.test <- t.test(cont.white$clintonid, treat.white$clintonid) #Compare whites in control versus treatment
print(white.test) #p=.006


#Compare blacks in treatment and control
cont.black.mean <- mean(cont.black$clintonid, na.rm = TRUE) #Get means and standard errors for blacks in control
cont.black.se <- sd(cont.black$clintonid, na.rm = TRUE)/sqrt(length(cont.black$clintonid))
print(cont.black.mean) #0.135
print(cont.black.se) #0.017

treat.black.mean <- mean(treat.black$clintonid, na.rm = TRUE) #Get means and standard errors for blacks in treatment
treat.black.se <- sd(treat.black$clintonid, na.rm = TRUE)/sqrt(length(treat.black$clintonid))
print(treat.black.mean) #0.212
print(treat.black.se) #0.020

black.test <- t.test(cont.black$clintonid, treat.black$clintonid) #Compare blacks in control versus treatment
print(black.test) #p=.004



#Compare white men in treatment and control
cont.white.men.mean <- mean(cont.white.men$clintonid, na.rm = TRUE) #Get means and standard errors for white men in control
cont.white.men.se <- sd(cont.white.men$clintonid, na.rm = TRUE)/sqrt(length(cont.white.men$clintonid))
print(cont.white.men.mean) #0.276
print(cont.white.men.se) #0.031

treat.white.men.mean <- mean(treat.white.men$clintonid, na.rm = TRUE) #Get means and standard errors for white men in treatment
treat.white.men.se <- sd(treat.white.men$clintonid, na.rm = TRUE)/sqrt(length(treat.white.men$clintonid))
print(treat.white.men.mean) #0.288
print(treat.white.men.se) #0.031

white.men.test <- t.test(cont.white.men$clintonid, treat.white.men$clintonid) #Compare white men in control versus treatment
print(white.men.test) #p=.781



#Compare white women in treatment and control
cont.white.women.mean <- mean(cont.white.women$clintonid, na.rm = TRUE) #Get means and standard errors for white women in control
cont.white.women.se <- sd(cont.white.women$clintonid, na.rm = TRUE)/sqrt(length(cont.white.women$clintonid))
print(cont.white.women.mean) #0.114
print(cont.white.women.se) #0.022

treat.white.women.mean <- mean(treat.white.women$clintonid, na.rm = TRUE) #Get means and standard errors for white women in control
treat.white.women.se <- sd(treat.white.women$clintonid, na.rm = TRUE)/sqrt(length(treat.white.women$clintonid))
print(treat.white.women.mean) #0.262
print(treat.white.women.se) #0.030

white.women.test <- t.test(cont.white.women$clintonid, treat.white.women$clintonid) #Compare white women in control versus treatment
print(white.women.test) # p=.000



#Compare Black men in treatment and control
cont.black.men.mean <- mean(cont.black.men$clintonid, na.rm = TRUE) #Get means and standard errors for Black men in control
cont.black.men.se <- sd(cont.black.men$clintonid, na.rm = TRUE)/sqrt(length(cont.black.men$clintonid))
print(cont.black.men.mean) #0.159
print(cont.black.men.se) #0.025

treat.black.men.mean <- mean(treat.black.men$clintonid, na.rm = TRUE) #Get means and standard errors for Black men in control
treat.black.men.se <- sd(treat.black.men$clintonid, na.rm = TRUE)/sqrt(length(treat.black.men$clintonid))
print(treat.black.men.mean) #0.257
print(treat.black.men.se) #0.031

black.men.test <- t.test(cont.black.men$clintonid, treat.black.men$clintonid) #Compare Black men in control versus treatment
print(black.men.test) #p=.015



#Compare Black women in treatment and control
cont.black.women.mean <- mean(cont.black.women$clintonid, na.rm = TRUE) #Get means and standard errors for Black women in control
cont.black.women.se <- sd(cont.black.women$clintonid, na.rm = TRUE)/sqrt(length(cont.black.women$clintonid))
print(cont.black.women.mean) #0.111
print(cont.black.women.se) #0.022

treat.black.women.mean <- mean(treat.black.women$clintonid, na.rm = TRUE) #Get means and standard errors for Black women in control
treat.black.women.se <- sd(treat.black.women$clintonid, na.rm = TRUE)/sqrt(length(treat.black.women$clintonid))
print(treat.black.women.mean) #0.171
print(treat.black.women.se) #0.026

black.women.test <- t.test(cont.black.women$clintonid, treat.black.women$clintonid) #Compare Black women in control versus treatment
print(black.women.test) #p=.080


#Combine means and standard errors into vectors, calculate confidence intervals
means <- c(cont.white.mean, treat.white.mean, cont.black.mean, treat.black.mean, cont.white.men.mean, treat.white.men.mean, cont.white.women.mean, treat.white.women.mean, cont.black.men.mean, treat.black.men.mean, cont.black.women.mean, treat.black.women.mean) * 100
ci.lower <- means - 1.96 * c(cont.white.se, treat.white.se, cont.black.se, treat.black.se, cont.white.men.se, treat.white.men.se, cont.white.women.se, treat.white.women.se, cont.black.men.se, treat.black.men.se, cont.black.women.se, treat.black.women.se) * 100
ci.upper <- means + 1.96 * c(cont.white.se, treat.white.se, cont.black.se, treat.black.se, cont.white.men.se, treat.white.men.se, cont.white.women.se, treat.white.women.se, cont.black.men.se, treat.black.men.se, cont.black.women.se, treat.black.women.se) * 100


#Graph means and and CI's
#Set dimensions to 6.5x4.75
par(mar = c(6, 3.5, 1, .25))
par(family = "serif")
x <- barplot(means, col = c("grey80", "grey40"), space = c(0, 0, 1, 0, 4, 0, 1, 0, 4, 0, 1, 0), ylim = c(0, 42), cex.axis = .9)
lines(c(-1, 23.5), c(0,0), xpd = NA)
text(c(1, 4, 10, 13, 19, 22), -2.75, c("All\nWhites", "All\nBlacks", "White\nMen", "White\nWomen", "Black\nMen", "Black\nWomen"), cex = .95, xpd = NA)
text(-3.20, 20, "Percentage", cex = .9, srt = 90, xpd = NA)
legend(-.35, -6.6, fill = c("grey80", "grey40"), c("Without Identity Politics Loss Narrative  ", "With Identity Politics Loss Narrative"), cex = .95, y.intersp = 1.1,  xpd = NA)
lines(c(x[1], x[1]), c(ci.lower[1], ci.upper[1]))
lines(c(x[1]-.05, x[1]+.05), c(ci.lower[1], ci.lower[1]))
lines(c(x[1]-.05, x[1]+.05), c(ci.upper[1], ci.upper[1]))
lines(c(x[2], x[2]), c(ci.lower[2], ci.upper[2]))
lines(c(x[2]-.05, x[2]+.05), c(ci.lower[2], ci.lower[2]))
lines(c(x[2]-.05, x[2]+.05), c(ci.upper[2], ci.upper[2]))
lines(c(x[3], x[3]), c(ci.lower[3], ci.upper[3]))
lines(c(x[3]-.05, x[3]+.05), c(ci.lower[3], ci.lower[3]))
lines(c(x[3]-.05, x[3]+.05), c(ci.upper[3], ci.upper[3]))
lines(c(x[4], x[4]), c(ci.lower[4], ci.upper[4]))
lines(c(x[4]-.05, x[4]+.05), c(ci.lower[4], ci.lower[4]))
lines(c(x[4]-.05, x[4]+.05), c(ci.upper[4], ci.upper[4]))
lines(c(x[5], x[5]), c(ci.lower[5], ci.upper[5]))
lines(c(x[5]-.05, x[5]+.05), c(ci.lower[5], ci.lower[5]))
lines(c(x[5]-.05, x[5]+.05), c(ci.upper[5], ci.upper[5]))
lines(c(x[6], x[6]), c(ci.lower[6], ci.upper[6]))
lines(c(x[6]-.05, x[6]+.05), c(ci.lower[6], ci.lower[6]))
lines(c(x[6]-.05, x[6]+.05), c(ci.upper[6], ci.upper[6]))
lines(c(x[7], x[7]), c(ci.lower[7], ci.upper[7]))
lines(c(x[7]-.05, x[7]+.05), c(ci.lower[7], ci.lower[7]))
lines(c(x[7]-.05, x[7]+.05), c(ci.upper[7], ci.upper[7]))
lines(c(x[8], x[8]), c(ci.lower[8], ci.upper[8]))
lines(c(x[8]-.05, x[8]+.05), c(ci.lower[8], ci.lower[8]))
lines(c(x[8]-.05, x[8]+.05), c(ci.upper[8], ci.upper[8]))
lines(c(x[9], x[9]), c(ci.lower[9], ci.upper[9]))
lines(c(x[9]-.05, x[9]+.05), c(ci.lower[9], ci.lower[9]))
lines(c(x[9]-.05, x[9]+.05), c(ci.upper[9], ci.upper[9]))
lines(c(x[10], x[10]), c(ci.lower[10], ci.upper[10]))
lines(c(x[10]-.05, x[10]+.05), c(ci.lower[10], ci.lower[10]))
lines(c(x[10]-.05, x[10]+.05), c(ci.upper[10], ci.upper[10]))
lines(c(x[11], x[11]), c(ci.lower[11], ci.upper[11]))
lines(c(x[11]-.05, x[11]+.05), c(ci.lower[11], ci.lower[11]))
lines(c(x[11]-.05, x[11]+.05), c(ci.upper[11], ci.upper[11]))
lines(c(x[12], x[12]), c(ci.lower[12], ci.upper[12]))
lines(c(x[12]-.05, x[12]+.05), c(ci.lower[12], ci.lower[12]))
lines(c(x[12]-.05, x[12]+.05), c(ci.upper[12], ci.upper[12]))
text(1, ci.upper[6] + 5, c("+8.0%\n p < .01"), cex = .85)
text(4, ci.upper[6] + 5, c("+7.7%\n p < .01"), cex = .85)
text(10, ci.upper[6] + 5, c("+1.2%\n p = .78"), cex = .85)
text(13, ci.upper[6] + 5, c("+14.7%\n p < .001"), cex = .85)
text(19, ci.upper[6] + 5, c("+9.9%\n p < .05"), cex = .85)
text(22, ci.upper[6] + 5, c("+6.0%\n p < .10"), cex = .85)
lines(c(x[1], x[1]), c(ci.upper[1] + 2, 37))
lines(c(x[2], x[2]), c(ci.upper[2] + 1, 37))
lines(c(x[1], x[2]), c(37, 37))
lines(c(x[3], x[3]), c(ci.upper[3] + 2, 37))
lines(c(x[4], x[4]), c(ci.upper[4] + 1, 37))
lines(c(x[3], x[4]), c(37, 37))
lines(c(x[5], x[5]), c(ci.upper[5] + 1, 37))
lines(c(x[6], x[6]), c(ci.upper[6] + 1, 37))
lines(c(x[5], x[6]), c(37, 37))
lines(c(x[7], x[7]), c(ci.upper[7] + 2, 37))
lines(c(x[8], x[8]), c(ci.upper[8] + 1, 37))
lines(c(x[7], x[8]), c(37, 37))
lines(c(x[9], x[9]), c(ci.upper[9] + 2, 37))
lines(c(x[10], x[10]), c(ci.upper[10] + 1, 37))
lines(c(x[9], x[10]), c(37, 37))
lines(c(x[11], x[11]), c(ci.upper[11] + 2, 37))
lines(c(x[12], x[12]), c(ci.upper[12] + 2, 37))
lines(c(x[11], x[12]), c(37, 37))



#####ORGANIZE DATA FOR REGRESSIONS#####

#Organize data for whites
smallstackedw <- subset(stacked.w, select = c(data.id, candvote, treatment, candgender, candrace, candideol, candpolicy, female)) #Create a subset that includes only the variables needed for the regressions
head(smallstackedw)
nrow(smallstackedw)

smallstackedw$candrace <- as.factor(smallstackedw$candrace) #Create factor variables for candidate traits - white, moderate, good jobs, male as the baseline
smallstackedw$candpolicy <- as.factor(smallstackedw$candpolicy)
smallstackedw$candideol <- as.factor(smallstackedw$candideol)

women.white <- subset(smallstackedw, female == 1) #Divide data between white men and women
men.white <- subset(smallstackedw, female == 0)


#Organize data for Blacks
smallstackedb <- subset(stacked.b, select = c(data.id, candvote, treatment, candgender, candrace, candideol, candpolicy,female)) #Create a subset that includes only the variables needed for the regressions
head(smallstackedb)
nrow(smallstackedb)

smallstackedb$candrace <- as.factor(smallstackedb$candrace) #Create factor variables candidate traits - white, moderate, good jobs, male as the baseline
smallstackedb$candpolicy <- as.factor(smallstackedb$candpolicy)
smallstackedb$candideol <- as.factor(smallstackedb$candideol)

women.black <- subset(smallstackedb, female == 1) #Divide data between Black men and women
men.black <- subset(smallstackedb, female == 0) #Divide data between Black men and women


#Create subsets by race and gender
women.black.control <- subset(women.black, treatment == 0)
women.black.treat <- subset(women.black, treatment == 1)
men.black.control <- subset(men.black, treatment == 0)
men.black.treat <- subset(men.black, treatment == 1)

women.white.control <- subset(women.white, treatment == 0)
women.white.treat <- subset(women.white, treatment == 1)
men.white.control <- subset(men.white, treatment == 0)
men.white.treat <- subset(men.white, treatment == 1)


#####REGRESSIONS FOR CANDIDATE VOTE CHOICE BY RACE AND GENDER (TABLE 2)#####

#Preliminary Regressions
wmc <- zelig(candvote ~ candgender + candrace + candideol + candpolicy, model = "logit.gee", id = "data.id", data = men.white.control, cite = FALSE) #White men in control
summary(wmc) #Ignore any error messages

wmt <- zelig(candvote ~ candgender + candrace + candideol + candpolicy, model = "logit.gee", id = "data.id", data = men.white.treat, cite = FALSE) #White men in treatment
summary(wmt)

wwc <- zelig(candvote ~ candgender + candrace + candideol + candpolicy, model = "logit.gee", id = "data.id", data = women.white.control, cite = FALSE) #White women in control
summary(wwc)

wwt <- zelig(candvote ~ candgender + candrace + candideol + candpolicy, model = "logit.gee", id = "data.id", data = women.white.treat, cite = FALSE) #White women in treatment
summary(wwt)


bmc <- zelig(candvote ~ candgender + candrace + candideol + candpolicy, model = "logit.gee", id = "data.id", data = men.black.control, cite = FALSE) #Black men in control
summary(bmc)

bmt <- zelig(candvote ~ candgender + candrace + candideol + candpolicy, model = "logit.gee", id = "data.id", data = men.black.treat, cite = FALSE) #Black men in treatment
summary(bmt)

bwc <- zelig(candvote ~ candgender + candrace + candideol + candpolicy, model = "logit.gee", id = "data.id", data = women.black.control, cite = FALSE) #Black women in control
summary(bwc)

bwt <- zelig(candvote ~ candgender + candrace + candideol + candpolicy, model = "logit.gee", id = "data.id", data = women.black.treat, cite = FALSE) #Black women in treatment
summary(bwt)


#Main Regressions
logit.white.men <- zelig(candvote ~ candgender*treatment + candrace*treatment + candideol*treatment + candpolicy*treatment, model = "logit.gee", id = "data.id", data = men.white, cite = FALSE) #White men
summary(logit.white.men)
nrow(na.omit(men.white))


logit.white.women <- zelig(candvote ~ candgender*treatment + candrace*treatment + candideol*treatment + candpolicy*treatment, model = "logit.gee", id = "data.id", data = women.white, cite = FALSE) #White women
summary(logit.white.women)
nrow(na.omit(women.white))


logit.black.men <- zelig(candvote ~ candgender*treatment + candrace*treatment + candideol*treatment + candpolicy*treatment, model = "logit.gee", id = "data.id", data = men.black, cite = FALSE)
summary(logit.black.men)
nrow(na.omit(men.black))


logit.black.women <- zelig(candvote ~ candgender*treatment + candrace*treatment + candideol*treatment + candpolicy*treatment, model = "logit.gee", id = "data.id", data = women.black, cite = FALSE)
summary(logit.black.women)
nrow(na.omit(women.black))


#####REGRESSIONS FOR MARGINAL PROBABILITY OF VOTE CHOICE (FIGURE 3)#####

#m1 is white men
#m2 is white women
#m3 is black men
#m4 is black women

#Get predicted probabilities for white men 
m1 <- zelig(candvote ~ candgender*treatment + candrace*treatment + candideol*treatment + candpolicy*treatment, model = "logit.gee", id = "data.id", data = men.white, cite = FALSE) #White men
summary(m1) 


#White men, gender probabilities
predprob1 <- setx(m1, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1)
predprob2 <- setx(m1, treatment = 0, candgender = 1, candrace = 1, candideol = 1, candpolicy = 1)
predprob3 <- setx(m1, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1)
predprob4 <- setx(m1, treatment = 1, candgender = 1, candrace = 1, candideol = 1, candpolicy = 1)

gender.control.whitemen <- sim(m1, x = predprob1, x1 = predprob2, num = 100000)
summary(gender.control.whitemen)

gender.treatment.whitemen <- sim(m1, x = predprob3, x1 = predprob4, num = 100000)
summary(gender.treatment.whitemen)

genderdiff.whitemen <- data.frame(matrix(ncol = 0, nrow = 2))
genderdiff.whitemen$point <- c(.040, .059) #Point estimates, may differ slightly by simulation
genderdiff.whitemen$bar <- genderdiff.whitemen$point - c(.008, .027) #Calculate length of lower CI
genderdiff.whitemen$label <- c("FemaleControl", "FemaleTreatment")
print(genderdiff.whitemen)

genderdid.sim <- (gender.treatment.whitemen$get_qi(qi = "fd", xvalue = "x1") -  gender.control.whitemen$get_qi(qi = "fd", xvalue = "x1")) #Calculate first difference bewteen treatment and control
summary(genderdid.sim)
sd(genderdid.sim)

genderdid.whitemen <- data.frame(matrix(ncol = 0, nrow = 1))
genderdid.whitemen$point <- .019 #Take men from genderdid.sim
genderdid.whitemen$bar <-.024*1.96 #Take sd from genderdid.sim
genderdid.whitemen$label <- "Female"
print(genderdid.whitemen)


#White men, race
predprob1 <- setx(m1, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1) #Get predicted probabilities for overmen race
predprob2 <- setx(m1, treatment = 0, candgender = 0, candrace = 2, candideol = 1, candpolicy = 1)
predprob3 <- setx(m1, treatment = 0, candgender = 0, candrace = 3, candideol = 1, candpolicy = 1)
predprob4 <- setx(m1, treatment = 0, candgender = 0, candrace = 4, candideol = 1, candpolicy = 1)
predprob5 <- setx(m1, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1)
predprob6 <- setx(m1, treatment = 1, candgender = 0, candrace = 2, candideol = 1, candpolicy = 1)
predprob7 <- setx(m1, treatment = 1, candgender = 0, candrace = 3, candideol = 1, candpolicy = 1)
predprob8 <- setx(m1, treatment = 1, candgender = 0, candrace = 4, candideol = 1, candpolicy = 1)

raceasian.control.whitemen <- sim(m1, x = predprob1, x1 = predprob2, num = 100000)
summary(raceasian.control.whitemen)

raceasian.treatment.whitemen <- sim(m1, x = predprob5, x1 = predprob6, num = 100000)
summary(raceasian.treatment.whitemen)

racehispanic.control.whitemen <- sim(m1, x = predprob1, x1 = predprob3, num = 100000)
summary(racehispanic.control.whitemen)

racehispanic.treatment.whitemen <- sim(m1, x = predprob5, x1 = predprob7, num = 100000)
summary(racehispanic.treatment.whitemen)

raceblack.control.whitemen <- sim(m1, x = predprob1, x1 = predprob4, num = 100000)
summary(raceblack.control.whitemen)

raceblack.treatment.whitemen <- sim(m1, x = predprob5, x1 = predprob8, num = 100000)
summary(raceblack.treatment.whitemen)

racediff.whitemen <- data.frame(matrix(ncol = 0, nrow = 6))
racediff.whitemen$point <- c(-.041, -.040, -.075, -.059, -.032, -.014)
racediff.whitemen$bar <- racediff.whitemen$point - c(-.085, -.082, -.123, -.105, -.078, -.059)
racediff.whitemen$label <- c("AsianControl", "AsianTreatment", "HispanicControl", "HispanicTreatment", "BlackControl", "BlackTreatment")
print(racediff.whitemen)

asian.sim <- (raceasian.treatment.whitemen$get_qi(qi = "fd", xvalue = "x1")) -  (raceasian.control.whitemen$get_qi(qi = "fd", xvalue = "x1"))
summary(asian.sim)
sd(asian.sim)

hisp.sim <- (racehispanic.treatment.whitemen$get_qi(qi = "fd", xvalue = "x1")) -  (racehispanic.control.whitemen$get_qi(qi = "fd", xvalue = "x1"))
summary(hisp.sim)
sd(hisp.sim)

black.sim <- (raceblack.treatment.whitemen$get_qi(qi = "fd", xvalue = "x1")) -  (raceblack.control.whitemen$get_qi(qi = "fd", xvalue = "x1"))
summary(black.sim)
sd(black.sim)

racedid.whitemen <- data.frame(matrix(ncol = 0, nrow = 3))
racedid.whitemen$point <- c(.001, .016, .018)
racedid.whitemen$bar <- c(.031, .035, .033)*1.96
racedid.whitemen$label <- c("Asian", "Hispanic", "Black")
print(racedid.whitemen)


#White men, ideology
predprob1 <- setx(m1, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1) #Get predicted probabilities for overmen ideology
predprob2 <- setx(m1, treatment = 0, candgender = 0, candrace = 1, candideol = 2, candpolicy = 1)
predprob3 <- setx(m1, treatment = 0, candgender = 0, candrace = 1, candideol = 3, candpolicy = 1)
predprob4 <- setx(m1, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1)
predprob5 <- setx(m1, treatment = 1, candgender = 0, candrace = 1, candideol = 2, candpolicy = 1)
predprob6 <- setx(m1, treatment = 1, candgender = 0, candrace = 1, candideol = 3, candpolicy = 1)

ideologycon.control.whitemen <- sim(m1, x = predprob1, x1 = predprob2, num = 100000)
summary(ideologycon.control.whitemen)

ideologycon.treatment.whitemen <- sim(m1, x = predprob4, x1 = predprob5, num = 100000)
summary(ideologycon.treatment.whitemen)

ideologylib.control.whitemen <- sim(m1, x = predprob1, x1 = predprob3, num = 100000)
summary(ideologylib.control.whitemen)

ideologylib.treatment.whitemen <- sim(m1, x = predprob4, x1 = predprob6, num = 100000)
summary(ideologylib.treatment.whitemen)

ideologydiff.whitemen <- data.frame(matrix(ncol = 0, nrow = 4))
ideologydiff.whitemen$point <- c(.021, .025, -.046, -.024)
ideologydiff.whitemen$bar <- ideologydiff.whitemen$point - c(-.022, -.014, -.095, -.074)
ideologydiff.whitemen$label <- c("ConservativeControl", "ConservativeTreatment", "LiberalControl", "LiberalTreatment")
print(ideologydiff.whitemen)

con.sim <- (ideologycon.treatment.whitemen$get_qi(qi = "fd", xvalue = "x1")) -  (ideologycon.control.whitemen$get_qi(qi = "fd", xvalue = "x1"))
summary(con.sim)
sd(con.sim)

lib.sim <- (ideologylib.treatment.whitemen$get_qi(qi = "fd", xvalue = "x1")) -  (ideologylib.control.whitemen$get_qi(qi = "fd", xvalue = "x1"))
summary(lib.sim)
sd(lib.sim)

ideologydid.whitemen <- data.frame(matrix(ncol = 0, nrow = 2))
ideologydid.whitemen$point <- c(.005, .022)
ideologydid.whitemen$bar <- c(.0295, .036)*1.96
ideologydid.whitemen$label <- c("Conservative", "Liberal")
print(ideologydid.whitemen)


#White men, policy
predprob1 <- setx(m1, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1)
predprob2 <- setx(m1, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 2)
predprob3 <- setx(m1, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 3)
predprob4 <- setx(m1, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 4)
predprob5 <- setx(m1, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1)
predprob6 <- setx(m1, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 2)
predprob7 <- setx(m1, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 3)
predprob8 <- setx(m1, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 4)


policyreform.control.whitemen <- sim(m1, x = predprob1, x1 = predprob2, num = 100000)
summary(policyreform.control.whitemen)

policyreform.treatment.whitemen <- sim(m1, x = predprob5, x1 = predprob6, num = 100000)
summary(policyreform.treatment.whitemen)

policyfairjob.control.whitemen <- sim(m1, x = predprob1, x1 = predprob3, num = 100000)
summary(policyfairjob.control.whitemen)

policyfairjob.treatment.whitemen <- sim(m1, x = predprob5, x1 = predprob7, num = 100000)
summary(policyfairjob.treatment.whitemen)

policyfairprison.control.whitemen <- sim(m1, x = predprob1, x1 = predprob4, num = 100000)
summary(policyfairprison.control.whitemen)

policyfairprison.treatment.whitemen <- sim(m1, x = predprob5, x1 = predprob8, num = 100000)
summary(policyfairprison.treatment.whitemen)


policydiff.whitemen <- data.frame(matrix(ncol = 0, nrow = 6))
policydiff.whitemen$point <- c(-.107, -.085, -.150, -.143, -.167, -.147)
policydiff.whitemen$bar <- policydiff.whitemen$point - c(-.159, -.136, -.205, -.193, -.228, -.203)
policydiff.whitemen$label <- c("GoodReformControl", "GoodReformTreatment", "FairJobsControl", "FairJobsTreatment", "FairReformControl", "FairReformTreatment")
print(policydiff.whitemen)

goodreform.sim <- (policyreform.treatment.whitemen$get_qi(qi = "fd", xvalue = "x1")) -  (policyreform.control.whitemen$get_qi(qi = "fd", xvalue = "x1"))
summary(goodreform.sim)
sd(goodreform.sim)

fairjobs.sim <- (policyfairjob.treatment.whitemen$get_qi(qi = "fd", xvalue = "x1")) -  (policyfairjob.control.whitemen$get_qi(qi = "fd", xvalue = "x1"))
summary(fairjobs.sim)
sd(fairjobs.sim)

fairreform.sim <- (policyfairprison.treatment.whitemen$get_qi(qi = "fd", xvalue = "x1")) -  (policyfairprison.control.whitemen$get_qi(qi = "fd", xvalue = "x1"))
summary(fairreform.sim)
sd(fairreform.sim)

policydid.whitemen <- data.frame(matrix(ncol = 0, nrow = 3))
policydid.whitemen$point <- c(.024, .009, .022)
policydid.whitemen$bar <- c(.038, .0386, .042)*1.96
policydid.whitemen$label <- c("GoodReform", "FairJobs", "FairReform")
print(policydid.whitemen)


whitemen.predprob <- rbind(genderdiff.whitemen, racediff.whitemen, ideologydiff.whitemen, policydiff.whitemen)
print(whitemen.predprob)

whitemen.did <- rbind(genderdid.whitemen, racedid.whitemen, ideologydid.whitemen, policydid.whitemen)
print(whitemen.did)



#Get predicted probabilities for white women
m2 <- zelig(candvote ~ candgender*treatment + candrace*treatment + candideol*treatment + candpolicy*treatment, model = "logit.gee", id = "data.id", data = women.white, cite = FALSE)
summary(m2)

#White women, gender
predprob1 <- setx(m2, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1) #Get predicted probabilities for white women for gender
predprob2 <- setx(m2, treatment = 0, candgender = 1, candrace = 1, candideol = 1, candpolicy = 1)
predprob3 <- setx(m2, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1)
predprob4 <- setx(m2, treatment = 1, candgender = 1, candrace = 1, candideol = 1, candpolicy = 1)

gender.control.whitewomen <- sim(m2, x = predprob1, x1 = predprob2, num = 100000)
summary(gender.control.whitewomen)

gender.treatment.whitewomen <- sim(m2, x = predprob3, x1 = predprob4, num = 100000)
summary(gender.treatment.whitewomen)

genderdiff.whitewomen <- data.frame(matrix(ncol = 0, nrow = 2))
genderdiff.whitewomen$point <- c(.113, .061)
genderdiff.whitewomen$bar <- genderdiff.whitewomen$point - c(.079, .033)
genderdiff.whitewomen$label <- c("FemaleControl", "FemaleTreatment")
print(genderdiff.whitewomen)

genderdid.sim <- (gender.treatment.whitewomen$get_qi(qi = "fd", xvalue = "x1")) -  (gender.control.whitewomen$get_qi(qi = "fd", xvalue = "x1"))
summary(genderdid.sim)
sd(genderdid.sim)

genderdid.whitewomen <- data.frame(matrix(ncol = 0, nrow = 1))
genderdid.whitewomen$point <- -.052
genderdid.whitewomen$bar <- .023*1.96
genderdid.whitewomen$label <- "Female"
print(genderdid.whitewomen)

#White women, race
predprob1 <- setx(m2, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1) #Get predicted probabilities for overwomen race
predprob2 <- setx(m2, treatment = 0, candgender = 0, candrace = 2, candideol = 1, candpolicy = 1)
predprob3 <- setx(m2, treatment = 0, candgender = 0, candrace = 3, candideol = 1, candpolicy = 1)
predprob4 <- setx(m2, treatment = 0, candgender = 0, candrace = 4, candideol = 1, candpolicy = 1)
predprob5 <- setx(m2, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1)
predprob6 <- setx(m2, treatment = 1, candgender = 0, candrace = 2, candideol = 1, candpolicy = 1)
predprob7 <- setx(m2, treatment = 1, candgender = 0, candrace = 3, candideol = 1, candpolicy = 1)
predprob8 <- setx(m2, treatment = 1, candgender = 0, candrace = 4, candideol = 1, candpolicy = 1)

raceasian.control.whitewomen <- sim(m2, x = predprob1, x1 = predprob2, num = 100000)
summary(raceasian.control.whitewomen)

raceasian.treatment.whitewomen <- sim(m2, x = predprob5, x1 = predprob6, num = 100000)
summary(raceasian.treatment.whitewomen)

racehispanic.control.whitewomen <- sim(m2, x = predprob1, x1 = predprob3, num = 100000)
summary(racehispanic.control.whitewomen)

racehispanic.treatment.whitewomen <- sim(m2, x = predprob5, x1 = predprob7, num = 100000)
summary(racehispanic.treatment.whitewomen)

raceblack.control.whitewomen <- sim(m2, x = predprob1, x1 = predprob4, num = 100000)
summary(raceblack.control.whitewomen)

raceblack.treatment.whitewomen <- sim(m2, x = predprob5, x1 = predprob8, num = 100000)
summary(raceblack.treatment.whitewomen)


racediff.whitewomen <- data.frame(matrix(ncol = 0, nrow = 6))
racediff.whitewomen$point <- c(-.042, -.022, -.018, -.037, .056, .014)
racediff.whitewomen$bar <- racediff.whitewomen$point - c(-.089, -.064, -.061, -.078, .011, -.024)
racediff.whitewomen$label <- c("AsianControl", "AsianTreatment", "HispanicControl", "HispanicTreatment", "BlackControl", "BlackTreatment")
print(racediff.whitewomen)


asian.sim <- (raceasian.treatment.whitewomen$get_qi(qi = "fd", xvalue = "x1")) -  (raceasian.control.whitewomen$get_qi(qi = "fd", xvalue = "x1"))
summary(asian.sim)
sd(asian.sim)

hisp.sim <- (racehispanic.treatment.whitewomen$get_qi(qi = "fd", xvalue = "x1")) -  (racehispanic.control.whitewomen$get_qi(qi = "fd", xvalue = "x1"))
summary(hisp.sim)
sd(hisp.sim)

black.sim <- (raceblack.treatment.whitewomen$get_qi(qi = "fd", xvalue = "x1")) -  (raceblack.control.whitewomen$get_qi(qi = "fd", xvalue = "x1"))
summary(black.sim)
sd(black.sim)

racedid.whitewomen <- data.frame(matrix(ncol = 0, nrow = 3))
racedid.whitewomen$point <- c(.018, -.021, -.041)
racedid.whitewomen$bar <- c(.032, .031, .031)*1.96
racedid.whitewomen$label <- c("Asian", "Hispanic", "Black")
print(racedid.whitewomen)


#White women, ideology
predprob1 <- setx(m2, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1) #Get predicted probabilities for overwomen ideology
predprob2 <- setx(m2, treatment = 0, candgender = 0, candrace = 1, candideol = 2, candpolicy = 1)
predprob3 <- setx(m2, treatment = 0, candgender = 0, candrace = 1, candideol = 3, candpolicy = 1)
predprob4 <- setx(m2, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1)
predprob5 <- setx(m2, treatment = 1, candgender = 0, candrace = 1, candideol = 2, candpolicy = 1)
predprob6 <- setx(m2, treatment = 1, candgender = 0, candrace = 1, candideol = 3, candpolicy = 1)

ideologycon.control.whitewomen <- sim(m2, x = predprob1, x1 = predprob2, num = 100000)
summary(ideologycon.control.whitewomen)

ideologycon.treatment.whitewomen <- sim(m2, x = predprob4, x1 = predprob5, num = 100000)
summary(ideologycon.treatment.whitewomen)

ideologylib.control.whitewomen <- sim(m2, x = predprob1, x1 = predprob3, num = 100000)
summary(ideologylib.control.whitewomen)

ideologylib.treatment.whitewomen <- sim(m2, x = predprob4, x1 = predprob6, num = 100000)
summary(ideologylib.treatment.whitewomen)


ideologydiff.whitewomen <- data.frame(matrix(ncol = 0, nrow = 4))
ideologydiff.whitewomen$point <- c(.055, .088, -.009, .006)
ideologydiff.whitewomen$bar <- ideologydiff.whitewomen$point - c(.014, .049, -.059, -.040)
ideologydiff.whitewomen$label <- c("ConservativeControl", "ConservativeTreatment", "LiberalControl", "LiberalTreatment")
print(ideologydiff.whitewomen)

con.sim <- (ideologycon.treatment.whitewomen$get_qi(qi = "fd", xvalue = "x1")) -  (ideologycon.control.whitewomen$get_qi(qi = "fd", xvalue = "x1"))
summary(con.sim)
sd(con.sim)

lib.sim <- (ideologylib.treatment.whitewomen$get_qi(qi = "fd", xvalue = "x1")) -  (ideologylib.control.whitewomen$get_qi(qi = "fd", xvalue = "x1"))
summary(lib.sim)
sd(lib.sim)

ideologydid.whitewomen <- data.frame(matrix(ncol = 0, nrow = 2))
ideologydid.whitewomen$point <- c(.033, .022)
ideologydid.whitewomen$bar <- c(.029, .036)*1.96
ideologydid.whitewomen$label <- c("Conservative", "Liberal")
print(ideologydid.whitewomen)


#White women, policy
predprob1 <- setx(m2, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1)
predprob2 <- setx(m2, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 2)
predprob3 <- setx(m2, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 3)
predprob4 <- setx(m2, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 4)
predprob5 <- setx(m2, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1)
predprob6 <- setx(m2, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 2)
predprob7 <- setx(m2, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 3)
predprob8 <- setx(m2, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 4)


policyreform.control.whitewomen <- sim(m2, x = predprob1, x1 = predprob2, num = 100000)
summary(policyreform.control.whitewomen)

policyreform.treatment.whitewomen <- sim(m2, x = predprob5, x1 = predprob6, num = 100000)
summary(policyreform.treatment.whitewomen)

policyfairjob.control.whitewomen <- sim(m2, x = predprob1, x1 = predprob3, num = 100000)
summary(policyfairjob.control.whitewomen)

policyfairjob.treatment.whitewomen <- sim(m2, x = predprob5, x1 = predprob7, num = 100000)
summary(policyfairjob.treatment.whitewomen)

policyfairprison.control.whitewomen <- sim(m2, x = predprob1, x1 = predprob4, num = 100000)
summary(policyfairprison.control.whitewomen)

policyfairprison.treatment.whitewomen <- sim(m2, x = predprob5, x1 = predprob8, num = 100000)
summary(policyfairprison.treatment.whitewomen)


policydiff.whitewomen <- data.frame(matrix(ncol = 0, nrow = 6))
policydiff.whitewomen$point <- c(-.078, -.151, -.009, -.11, -.115, -.225)
policydiff.whitewomen$bar <- policydiff.whitewomen$point - c(-.142, -.206, -.073, -.167, -.186, -.288)
policydiff.whitewomen$label <- c("GoodReformControl", "GoodReformTreatment", "FairJobsControl", "FairJobsTreatment", "FairReformControl", "FairReformTreatment")
print(policydiff.whitewomen)


goodreform.sim <- (policyreform.treatment.whitewomen$get_qi(qi = "fd", xvalue = "x1")) -  (policyreform.control.whitewomen$get_qi(qi = "fd", xvalue = "x1"))
summary(goodreform.sim)
sd(goodreform.sim)

fairjobs.sim <- (policyfairjob.treatment.whitewomen$get_qi(qi = "fd", xvalue = "x1")) -  (policyfairjob.control.whitewomen$get_qi(qi = "fd", xvalue = "x1"))
summary(fairjobs.sim)
sd(fairjobs.sim)

fairreform.sim <- (policyfairprison.treatment.whitewomen$get_qi(qi = "fd", xvalue = "x1")) -  (policyfairprison.control.whitewomen$get_qi(qi = "fd", xvalue = "x1"))
summary(fairreform.sim)
sd(fairreform.sim)

policydid.whitewomen <- data.frame(matrix(ncol = 0, nrow = 3))
policydid.whitewomen$point <- c(-.079, -.105, -.111)
policydid.whitewomen$bar <- c(.043, .045, .048)*1.96
policydid.whitewomen$label <- c("GoodReform", "FairJobs", "FairReform")
print(policydid.whitewomen)


whitewomen.predprob <- rbind(genderdiff.whitewomen, racediff.whitewomen, ideologydiff.whitewomen, policydiff.whitewomen)
print(whitewomen.predprob)

whitewomen.did <- rbind(genderdid.whitewomen, racedid.whitewomen, ideologydid.whitewomen, policydid.whitewomen)
print(whitewomen.did)


#Get predicted probabilities for black men
m3 <- zelig(candvote ~ candgender*treatment + candrace*treatment + candideol*treatment + candpolicy*treatment, model = "logit.gee", id = "data.id", data = men.black, cite = FALSE)
summary(m3)

#Get predicted probabilities for men respondents

#Men respondents, gender
predprob1 <- setx(m3, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1) #Get predicted probabilities for overmen gender
predprob2 <- setx(m3, treatment = 0, candgender = 1, candrace = 1, candideol = 1, candpolicy = 1)
predprob3 <- setx(m3, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1)
predprob4 <- setx(m3, treatment = 1, candgender = 1, candrace = 1, candideol = 1, candpolicy = 1)

gender.control.blackmen <- sim(m3, x = predprob1, x1 = predprob2, num = 100000)
summary(gender.control.blackmen)

gender.treatment.blackmen <- sim(m3, x = predprob3, x1 = predprob4, num = 100000)
summary(gender.treatment.blackmen)

genderdiff.blackmen <- data.frame(matrix(ncol = 0, nrow = 2))
genderdiff.blackmen$point <- c(-.003, .001)
genderdiff.blackmen$bar <- genderdiff.blackmen$point - c(-.038, -.038)
genderdiff.blackmen$label <- c("FemaleControl", "FemaleTreatment")
print(genderdiff.blackmen)

genderdid.sim <- (gender.treatment.blackmen$get_qi(qi = "fd", xvalue = "x1")) -  (gender.control.blackmen$get_qi(qi = "fd", xvalue = "x1"))
summary(genderdid.sim)
sd(genderdid.sim)

genderdid.blackmen <- data.frame(matrix(ncol = 0, nrow = 1))
genderdid.blackmen$point <- .0038
genderdid.blackmen$bar <-.026 * 1.96
genderdid.blackmen$label <- "Female"
print(genderdid.blackmen)


#Black men, race
predprob1 <- setx(m3, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1) #Get predicted probabilities for overmen race
predprob2 <- setx(m3, treatment = 0, candgender = 0, candrace = 2, candideol = 1, candpolicy = 1)
predprob3 <- setx(m3, treatment = 0, candgender = 0, candrace = 3, candideol = 1, candpolicy = 1)
predprob4 <- setx(m3, treatment = 0, candgender = 0, candrace = 4, candideol = 1, candpolicy = 1)
predprob5 <- setx(m3, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1)
predprob6 <- setx(m3, treatment = 1, candgender = 0, candrace = 2, candideol = 1, candpolicy = 1)
predprob7 <- setx(m3, treatment = 1, candgender = 0, candrace = 3, candideol = 1, candpolicy = 1)
predprob8 <- setx(m3, treatment = 1, candgender = 0, candrace = 4, candideol = 1, candpolicy = 1)

raceasian.control.blackmen <- sim(m3, x = predprob1, x1 = predprob2, num = 100000)
summary(raceasian.control.blackmen)

raceasian.treatment.blackmen <- sim(m3, x = predprob5, x1 = predprob6, num = 100000)
summary(raceasian.treatment.blackmen)

racehispanic.control.blackmen <- sim(m3, x = predprob1, x1 = predprob3, num = 100000)
summary(racehispanic.control.blackmen)

racehispanic.treatment.blackmen <- sim(m3, x = predprob5, x1 = predprob7, num = 100000)
summary(racehispanic.treatment.blackmen)

raceblack.control.blackmen <- sim(m3, x = predprob1, x1 = predprob4, num = 100000)
summary(raceblack.control.blackmen)

raceblack.treatment.blackmen <- sim(m3, x = predprob5, x1 = predprob8, num = 100000)
summary(raceblack.treatment.blackmen)


racediff.blackmen <- data.frame(matrix(ncol = 0, nrow = 6))
racediff.blackmen$point <- c(.011, .058, .023, .089, .178, .201)
racediff.blackmen$bar <- racediff.blackmen$point - c(-.033, .008, -.023, .043, .135, .156)
racediff.blackmen$label <- c("AsianControl", "AsianTreatment", "HispanicControl", "HispanicTreatment", "BlackControl", "BlackTreatment")
print(racediff.blackmen)

asian.sim <- (raceasian.treatment.blackmen$get_qi(qi = "fd", xvalue = "x1")) -  (raceasian.control.blackmen$get_qi(qi = "fd", xvalue = "x1"))
summary(asian.sim)
sd(asian.sim)

hisp.sim <- (racehispanic.treatment.blackmen$get_qi(qi = "fd", xvalue = "x1")) -  (racehispanic.control.blackmen$get_qi(qi = "fd", xvalue = "x1"))
summary(hisp.sim)
sd(hisp.sim)

black.sim <- (raceblack.treatment.blackmen$get_qi(qi = "fd", xvalue = "x1")) -  (raceblack.control.blackmen$get_qi(qi = "fd", xvalue = "x1"))
summary(black.sim)
sd(black.sim)

racedid.blackmen <- data.frame(matrix(ncol = 0, nrow = 3))
racedid.blackmen$point <- c(.047, .067, .024)
racedid.blackmen$bar <- c(.034, .033, .032)*1.96
racedid.blackmen$label <- c("Asian", "Hispanic", "Black")
print(racedid.blackmen)


#Black men, ideology
predprob1 <- setx(m3, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1) #Get predicted probabilities for overmen ideology
predprob2 <- setx(m3, treatment = 0, candgender = 0, candrace = 1, candideol = 2, candpolicy = 1)
predprob3 <- setx(m3, treatment = 0, candgender = 0, candrace = 1, candideol = 3, candpolicy = 1)
predprob4 <- setx(m3, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1)
predprob5 <- setx(m3, treatment = 1, candgender = 0, candrace = 1, candideol = 2, candpolicy = 1)
predprob6 <- setx(m3, treatment = 1, candgender = 0, candrace = 1, candideol = 3, candpolicy = 1)

ideologycon.control.blackmen <- sim(m3, x = predprob1, x1 = predprob2, num = 100000)
summary(ideologycon.control.blackmen)

ideologycon.treatment.blackmen <- sim(m3, x = predprob4, x1 = predprob5, num = 100000)
summary(ideologycon.treatment.blackmen)

ideologylib.control.blackmen <- sim(m3, x = predprob1, x1 = predprob3, num = 100000)
summary(ideologylib.control.blackmen)

ideologylib.treatment.blackmen <- sim(m3, x = predprob4, x1 = predprob6, num = 100000)
summary(ideologylib.treatment.blackmen)


ideologydiff.blackmen <- data.frame(matrix(ncol = 0, nrow = 4))
ideologydiff.blackmen$point <- c(-.025, -.025, .012, -0.019)
ideologydiff.blackmen$bar <- ideologydiff.blackmen$point - c(-.061, -.065, -.032, -.056)
ideologydiff.blackmen$label <- c("ConservativeControl", "ConservativeTreatment", "LiberalControl", "LiberalTreatment")
print(ideologydiff.blackmen)

con.sim <- (ideologycon.treatment.blackmen$get_qi(qi = "fd", xvalue = "x1")) -  (ideologycon.control.blackmen$get_qi(qi = "fd", xvalue = "x1"))
summary(con.sim)
sd(con.sim)

lib.sim <- (ideologylib.treatment.blackmen$get_qi(qi = "fd", xvalue = "x1")) -  (ideologylib.control.blackmen$get_qi(qi = "fd", xvalue = "x1"))
summary(lib.sim)
sd(lib.sim)

ideologydid.blackmen <- data.frame(matrix(ncol = 0, nrow = 2))
ideologydid.blackmen$point <- c(-.0002, 0.031)
ideologydid.blackmen$bar <- c(.027, .029)*1.96
ideologydid.blackmen$label <- c("Conservative", "Liberal")
print(ideologydid.blackmen)


#Black men, policy
predprob1 <- setx(m3, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1)
predprob2 <- setx(m3, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 2)
predprob3 <- setx(m3, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 3)
predprob4 <- setx(m3, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 4)
predprob5 <- setx(m3, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1)
predprob6 <- setx(m3, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 2)
predprob7 <- setx(m3, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 3)
predprob8 <- setx(m3, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 4)


policyreform.control.blackmen <- sim(m3, x = predprob1, x1 = predprob2, num = 100000)
summary(policyreform.control.blackmen)

policyreform.treatment.blackmen <- sim(m3, x = predprob5, x1 = predprob6, num = 100000)
summary(policyreform.treatment.blackmen)

policyfairjob.control.blackmen <- sim(m3, x = predprob1, x1 = predprob3, num = 100000)
summary(policyfairjob.control.blackmen)

policyfairjob.treatment.blackmen <- sim(m3, x = predprob5, x1 = predprob7, num = 100000)
summary(policyfairjob.treatment.blackmen)

policyfairprison.control.blackmen <- sim(m3, x = predprob1, x1 = predprob4, num = 100000)
summary(policyfairprison.control.blackmen)

policyfairprison.treatment.blackmen <- sim(m3, x = predprob5, x1 = predprob8, num = 100000)
summary(policyfairprison.treatment.blackmen)


policydiff.blackmen <- data.frame(matrix(ncol = 0, nrow = 6))
policydiff.blackmen$point <- c(-.014, -.087, .037, -.030, .089, .014)
policydiff.blackmen$bar <- policydiff.blackmen$point - c(-.066, -.141, -.019, -.081, .033, -.038)
policydiff.blackmen$label <- c("GoodReformControl", "GoodReformTreatment", "FairJobsControl", "FairJobsTreatment", "FairReformControl", "FairReformTreatment")
print(policydiff.blackmen)

goodreform.sim <- (policyreform.treatment.blackmen$get_qi(qi = "fd", xvalue = "x1")) -  (policyreform.control.blackmen$get_qi(qi = "fd", xvalue = "x1"))
summary(goodreform.sim)
sd(goodreform.sim)

fairjobs.sim <- (policyfairjob.treatment.blackmen$get_qi(qi = "fd", xvalue = "x1")) -  (policyfairjob.control.blackmen$get_qi(qi = "fd", xvalue = "x1"))
summary(fairjobs.sim)
sd(fairjobs.sim)

fairreform.sim <- (policyfairprison.treatment.blackmen$get_qi(qi = "fd", xvalue = "x1")) -  (policyfairprison.control.blackmen$get_qi(qi = "fd", xvalue = "x1"))
summary(fairreform.sim)
sd(fairreform.sim)

policydid.blackmen <- data.frame(matrix(ncol = 0, nrow = 3))
policydid.blackmen$point <- c(-.073, -.067, -.075)
policydid.blackmen$bar <- c(.038, .038, .039)*1.96
policydid.blackmen$label <- c("GoodReform", "FairJobs", "FairReform")
print(policydid.blackmen)


blackmen.predprob <- rbind(genderdiff.blackmen, racediff.blackmen, ideologydiff.blackmen, policydiff.blackmen)
print(blackmen.predprob)

blackmen.did <- rbind(genderdid.blackmen, racedid.blackmen, ideologydid.blackmen, policydid.blackmen)
print(blackmen.did)



#Get predicted probabilities for black women

m4 <- zelig(candvote ~ candgender*treatment + candrace*treatment + candideol*treatment + candpolicy*treatment, model = "logit.gee", id = "data.id", data = women.black, cite = FALSE)
summary(m4)

#Black women, gender
predprob1 <- setx(m4, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1) #Get predicted probabilities for overwomen gender
predprob2 <- setx(m4, treatment = 0, candgender = 1, candrace = 1, candideol = 1, candpolicy = 1)
predprob3 <- setx(m4, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1)
predprob4 <- setx(m4, treatment = 1, candgender = 1, candrace = 1, candideol = 1, candpolicy = 1)

gender.control.blackwomen <- sim(m4, x = predprob1, x1 = predprob2, num = 100000)
summary(gender.control.blackwomen)

gender.treatment.blackwomen <- sim(m4, x = predprob3, x1 = predprob4, num = 100000)
summary(gender.treatment.blackwomen)

genderdiff.blackwomen <- data.frame(matrix(ncol = 0, nrow = 2))
genderdiff.blackwomen$point <- c(.038, .033)
genderdiff.blackwomen$bar <- genderdiff.blackwomen$point - c(.003, -.002)
genderdiff.blackwomen$label <- c("FemaleControl", "FemaleTreatment")
print(genderdiff.blackwomen)

genderdid.sim <- (gender.treatment.blackwomen$get_qi(qi = "fd", xvalue = "x1")) -  (gender.control.blackwomen$get_qi(qi = "fd", xvalue = "x1"))
summary(genderdid.sim)
sd(genderdid.sim)

genderdid.blackwomen <- data.frame(matrix(ncol = 0, nrow = 1))
genderdid.blackwomen$point <- -.006
genderdid.blackwomen$bar <- .025*1.96
genderdid.blackwomen$label <- "Female"
print(genderdid.blackwomen)

#Black women, race
predprob1 <- setx(m4, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1) #Get predicted probabilities for overwomen race
predprob2 <- setx(m4, treatment = 0, candgender = 0, candrace = 2, candideol = 1, candpolicy = 1)
predprob3 <- setx(m4, treatment = 0, candgender = 0, candrace = 3, candideol = 1, candpolicy = 1)
predprob4 <- setx(m4, treatment = 0, candgender = 0, candrace = 4, candideol = 1, candpolicy = 1)
predprob5 <- setx(m4, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1)
predprob6 <- setx(m4, treatment = 1, candgender = 0, candrace = 2, candideol = 1, candpolicy = 1)
predprob7 <- setx(m4, treatment = 1, candgender = 0, candrace = 3, candideol = 1, candpolicy = 1)
predprob8 <- setx(m4, treatment = 1, candgender = 0, candrace = 4, candideol = 1, candpolicy = 1)

raceasian.control.blackwomen <- sim(m4, x = predprob1, x1 = predprob2, num = 100000)
summary(raceasian.control.blackwomen)

raceasian.treatment.blackwomen <- sim(m4, x = predprob5, x1 = predprob6, num = 100000)
summary(raceasian.treatment.blackwomen)

racehispanic.control.blackwomen <- sim(m4, x = predprob1, x1 = predprob3, num = 100000)
summary(racehispanic.control.blackwomen)

racehispanic.treatment.blackwomen <- sim(m4, x = predprob5, x1 = predprob7, num = 100000)
summary(racehispanic.treatment.blackwomen)

raceblack.control.blackwomen <- sim(m4, x = predprob1, x1 = predprob4, num = 100000)
summary(raceblack.control.blackwomen)

raceblack.treatment.blackwomen <- sim(m4, x = predprob5, x1 = predprob8, num = 100000)
summary(raceblack.treatment.blackwomen)


racediff.blackwomen <- data.frame(matrix(ncol = 0, nrow = 6))
racediff.blackwomen$point <- c(.025, -.006, .073, .027, .205, .205)
racediff.blackwomen$bar <- racediff.blackwomen$point - c(-.022, -.050, .024, -.015, .157, .158)
racediff.blackwomen$label <- c("AsianControl", "AsianTreatment", "HispanicControl", "HispanicTreatment", "BlackControl", "BlackTreatment")
print(racediff.blackwomen)


asian.sim <- (raceasian.treatment.blackwomen$get_qi(qi = "fd", xvalue = "x1")) -  (raceasian.control.blackwomen$get_qi(qi = "fd", xvalue = "x1"))
summary(asian.sim)
sd(asian.sim)

hisp.sim <- (racehispanic.treatment.blackwomen$get_qi(qi = "fd", xvalue = "x1")) -  (racehispanic.control.blackmen$get_qi(qi = "fd", xvalue = "x1"))
summary(hisp.sim)
sd(hisp.sim)

black.sim <- (raceblack.treatment.blackwomen$get_qi(qi = "fd", xvalue = "x1")) -  (raceblack.control.blackmen$get_qi(qi = "fd", xvalue = "x1"))
summary(black.sim)
sd(black.sim)

racedid.blackwomen <- data.frame(matrix(ncol = 0, nrow = 3))
racedid.blackwomen$point <- c(-.031, .005, .027)
racedid.blackwomen$bar <- c(.033, .031, .033)*1.96
racedid.blackwomen$label <- c("Asian", "Hispanic", "Black")
print(racedid.blackwomen)


#Black women, ideology
predprob1 <- setx(m4, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1) #Get predicted probabilities for overwomen ideology
predprob2 <- setx(m4, treatment = 0, candgender = 0, candrace = 1, candideol = 2, candpolicy = 1)
predprob3 <- setx(m4, treatment = 0, candgender = 0, candrace = 1, candideol = 3, candpolicy = 1)
predprob4 <- setx(m4, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1)
predprob5 <- setx(m4, treatment = 1, candgender = 0, candrace = 1, candideol = 2, candpolicy = 1)
predprob6 <- setx(m4, treatment = 1, candgender = 0, candrace = 1, candideol = 3, candpolicy = 1)

ideologycon.control.blackwomen <- sim(m4, x = predprob1, x1 = predprob2, num = 100000)
summary(ideologycon.control.blackwomen)

ideologycon.treatment.blackwomen <- sim(m4, x = predprob4, x1 = predprob5, num = 100000)
summary(ideologycon.treatment.blackwomen)

ideologylib.control.blackwomen <- sim(m4, x = predprob1, x1 = predprob3, num = 100000)
summary(ideologylib.control.blackwomen)

ideologylib.treatment.blackwomen <- sim(m4, x = predprob4, x1 = predprob6, num = 100000)
summary(ideologylib.treatment.blackwomen)


ideologydiff.blackwomen <- data.frame(matrix(ncol = 0, nrow = 4))
ideologydiff.blackwomen$point <- c(-.070, -.068, -.021, -.019)
ideologydiff.blackwomen$bar <- ideologydiff.blackwomen$point - c(-.107, -.105, -.062, -.056)
ideologydiff.blackwomen$label <- c("ConservativeControl", "ConservativeTreatment", "LiberalControl", "LiberalTreatment")
print(ideologydiff.blackwomen)

con.sim <- (ideologycon.treatment.blackwomen$get_qi(qi = "fd", xvalue = "x1")) -  (ideologycon.control.blackwomen$get_qi(qi = "fd", xvalue = "x1"))
summary(con.sim)
sd(con.sim)

lib.sim <- (ideologylib.treatment.blackwomen$get_qi(qi = "fd", xvalue = "x1")) -  (ideologylib.control.blackwomen$get_qi(qi = "fd", xvalue = "x1"))
summary(lib.sim)
sd(lib.sim)

ideologydid.blackwomen <- data.frame(matrix(ncol = 0, nrow = 2))
ideologydid.blackwomen$point <- c(.002, 0.003)
ideologydid.blackwomen$bar <- c(.027, .029)*1.96
ideologydid.blackwomen$label <- c("Conservative", "Liberal")
print(ideologydid.blackwomen)


#Black women, policy
predprob1 <- setx(m4, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1)
predprob2 <- setx(m4, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 2)
predprob3 <- setx(m4, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 3)
predprob4 <- setx(m4, treatment = 0, candgender = 0, candrace = 1, candideol = 1, candpolicy = 4)
predprob5 <- setx(m4, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 1)
predprob6 <- setx(m4, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 2)
predprob7 <- setx(m4, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 3)
predprob8 <- setx(m4, treatment = 1, candgender = 0, candrace = 1, candideol = 1, candpolicy = 4)


policyreform.control.blackwomen <- sim(m4, x = predprob1, x1 = predprob2, num = 100000)
summary(policyreform.control.blackwomen)

policyreform.treatment.blackwomen <- sim(m4, x = predprob5, x1 = predprob6, num = 100000)
summary(policyreform.treatment.blackwomen)

policyfairjob.control.blackwomen <- sim(m4, x = predprob1, x1 = predprob3, num = 100000)
summary(policyfairjob.control.blackwomen)

policyfairjob.treatment.blackwomen <- sim(m4, x = predprob5, x1 = predprob7, num = 100000)
summary(policyfairjob.treatment.blackwomen)

policyfairprison.control.blackwomen <- sim(m4, x = predprob1, x1 = predprob4, num = 100000)
summary(policyfairprison.control.blackwomen)

policyfairprison.treatment.blackwomen <- sim(m4, x = predprob5, x1 = predprob8, num = 100000)
summary(policyfairprison.treatment.blackwomen)


policydiff.blackwomen <- data.frame(matrix(ncol = 0, nrow = 6))
policydiff.blackwomen$point <- c(-.0036, -.007, .079, .031, .082, .075)
policydiff.blackwomen$bar <- policydiff.blackwomen$point - c(-.057, -.061, .026, -.024, .028, .019)
policydiff.blackwomen$label <- c("GoodReformControl", "GoodReformTreatment", "FairJobsControl", "FairJobsTreatment", "FairReformControl", "FairReformTreatment")
print(policydiff.blackwomen)


goodreform.sim <- (policyreform.treatment.blackwomen$get_qi(qi = "fd", xvalue = "x1")) -  (policyreform.control.blackwomen$get_qi(qi = "fd", xvalue = "x1"))
summary(goodreform.sim)
sd(goodreform.sim)

fairjobs.sim <- (policyfairjob.treatment.blackwomen$get_qi(qi = "fd", xvalue = "x1")) -  (policyfairjob.control.blackwomen$get_qi(qi = "fd", xvalue = "x1"))
summary(fairjobs.sim)
sd(fairjobs.sim)

fairreform.sim <- (policyfairprison.treatment.blackwomen$get_qi(qi = "fd", xvalue = "x1")) -  (policyfairprison.control.blackwomen$get_qi(qi = "fd", xvalue = "x1"))
summary(fairreform.sim)
sd(fairreform.sim)


policydid.blackwomen <- data.frame(matrix(ncol = 0, nrow = 3))
policydid.blackwomen$point <- c(-.004, -.047, -.007)
policydid.blackwomen$bar <- c(.039, .039, .040)*1.96
policydid.blackwomen$label <- c("GoodReform", "FairJobs", "FairReform")
print(policydid.blackwomen)


blackwomen.predprob <- rbind(genderdiff.blackwomen, racediff.blackwomen, ideologydiff.blackwomen, policydiff.blackwomen)
print(blackwomen.predprob)

blackwomen.did <- rbind(genderdid.blackwomen, racedid.blackwomen, ideologydid.blackwomen, policydid.blackwomen)
print(blackwomen.did)


#Create graph

fullyvalue <- rev(c(1, 1.5, 2.2, 4.2, 4.7, 5.4, 7.4, 7.9, 8.6,9.6, 12.6, 13.1, 13.8,15.8, 16.3, 17.0, 18.0, 21.0, 21.5, 22.2, 24.2, 24.7, 25.4, 27.4, 27.9, 28.6, 29.6, 32.6, 33.1, 33.8, 34.8))
pointyvalue <- rev(c(1, 1.5, 4.2, 4.7, 7.4, 7.9, 12.6, 13.1, 15.8, 16.3, 21.0, 21.5, 24.2, 24.7, 27.4, 27.9, 32.6, 33.1))
groupyvalue <- c(2.2, 5.4, 8.6, 13.8, 17.0, 22.2, 25.4, 28.6, 33.8)
groupyvalue.rev <- rev(c(2.2, 5.4, 8.6, 13.8, 17.0, 22.2, 25.4, 28.6, 33.8))
labyvalue <- c(9.6, 18.0, 29.6, 34.8)

whitemen.predprob$yvalue <- pointyvalue
whitemen.predprob$upper <- whitemen.predprob$point + whitemen.predprob$bar
whitemen.predprob$lower <- whitemen.predprob$point - whitemen.predprob$bar
print(whitemen.predprob)

whitewomen.predprob$yvalue <- pointyvalue
whitewomen.predprob$upper <- whitewomen.predprob$point + whitewomen.predprob$bar
whitewomen.predprob$lower <- whitewomen.predprob$point - whitewomen.predprob$bar
print(whitewomen.predprob)

blackmen.predprob$yvalue <- pointyvalue
blackmen.predprob$upper <- blackmen.predprob$point + blackmen.predprob$bar
blackmen.predprob$lower <- blackmen.predprob$point - blackmen.predprob$bar
print(blackmen.predprob)

blackwomen.predprob$yvalue <- pointyvalue
blackwomen.predprob$upper <- blackwomen.predprob$point + blackwomen.predprob$bar
blackwomen.predprob$lower <- blackwomen.predprob$point - blackwomen.predprob$bar
print(blackwomen.predprob)


par(mar=c(3, .45, 2, .45), oma = c(5, 11, 0, 0))
par(family = "serif")
par(mfrow = c(1,4))

plotCI(whitemen.predprob$point, whitemen.predprob$yvalue, xlab = "", ui = whitemen.predprob$upper, li = whitemen.predprob$lower, sfrac = .002, err = "x", lwd = 1, gap = 0, pch = c(21, 16), yaxt = "n", ylab = "", main = "White Men", xlim = c(-.31, .31), ylim = c(.5, 35.3), cex.axis = .8, cex.main = 1)
text(-.33, pointyvalue, pos = 2, c("Control", "Treatment", "Control", "Treatment","Control", "Treatment","Control", "Treatment","Control", "Treatment","Control", "Treatment") , srt = 360, xpd = NA)
text(-.33, groupyvalue, c("Fair Sentencing", "Fair Employment", "Justice Reform", "Very Liberal", "Slightly Liberal", "African American", "Hispanic", "Asian American", "Female"), pos = 2, font = 3, xpd = NA)
text(-.33, labyvalue, c("Policy (Good Jobs Baseline)", "Ideology (Moderate Baseline)", "Race (White Basline)", "Gender (Male Baseline)"), pos = 2, xpd = NA, font = 2)
abline(v = 0, lty = 2)


plotCI(whitewomen.predprob$point, whitewomen.predprob$yvalue, xlab = "", ui = whitewomen.predprob$upper, li = whitewomen.predprob$lower, sfrac = .002, err = "x", lwd = 1, gap = 0, pch = c(21, 8, 21, 16, 21, 16, 21, 16, 21, 16, 21, 16, 21, 8, 21, 8, 21, 8), yaxt = "n", ylab = "", main = "White Women", xlim = c(-.31, .31), ylim = c(.5, 35.3), cex.axis = .8, cex.main = 1)
abline(v = 0, lty = 2)

lines(c(whitewomen.predprob$upper[1] + .05, whitewomen.predprob$upper[1] + .05), c(whitewomen.predprob$yvalue[1], whitewomen.predprob$yvalue[2]))
lines(c(whitewomen.predprob$upper[1] + .05, whitewomen.predprob$upper[1] + .02), c(whitewomen.predprob$yvalue[1], whitewomen.predprob$yvalue[1]))
lines(c(whitewomen.predprob$upper[1] + .05, whitewomen.predprob$upper[1] + .02), c(whitewomen.predprob$yvalue[2], whitewomen.predprob$yvalue[2]))
text(whitewomen.predprob$upper[1] + .12, 32.9, "-.05", cex = .9)

lines(c(whitewomen.predprob$upper[13] + .05, whitewomen.predprob$upper[13] + .05), c(whitewomen.predprob$yvalue[13], whitewomen.predprob$yvalue[14]))
lines(c(whitewomen.predprob$upper[13] + .05, whitewomen.predprob$upper[13] + .02), c(whitewomen.predprob$yvalue[13], whitewomen.predprob$yvalue[13]))
lines(c(whitewomen.predprob$upper[13] + .05, whitewomen.predprob$upper[13] + .02), c(whitewomen.predprob$yvalue[14], whitewomen.predprob$yvalue[14]))
text(whitewomen.predprob$upper[13] + .12, 7.7, "-.07", cex = .9)

lines(c(whitewomen.predprob$upper[15] + .05, whitewomen.predprob$upper[15] + .05), c(whitewomen.predprob$yvalue[15], whitewomen.predprob$yvalue[16]))
lines(c(whitewomen.predprob$upper[15] + .05, whitewomen.predprob$upper[15] + .02), c(whitewomen.predprob$yvalue[15], whitewomen.predprob$yvalue[15]))
lines(c(whitewomen.predprob$upper[15] + .05, whitewomen.predprob$upper[15] + .02), c(whitewomen.predprob$yvalue[16], whitewomen.predprob$yvalue[16]))
text(whitewomen.predprob$upper[15] + .12, 4.5, "-.10", cex = .9)

lines(c(whitewomen.predprob$upper[17] + .09, whitewomen.predprob$upper[17] + .09), c(whitewomen.predprob$yvalue[17], whitewomen.predprob$yvalue[18]))
lines(c(whitewomen.predprob$upper[17] + .09, whitewomen.predprob$upper[17] + .06), c(whitewomen.predprob$yvalue[17], whitewomen.predprob$yvalue[17]))
lines(c(whitewomen.predprob$upper[17] + .09, whitewomen.predprob$upper[17] + .06), c(whitewomen.predprob$yvalue[18], whitewomen.predprob$yvalue[18]))
text(whitewomen.predprob$upper[17] + .16, 1.3, "-.11", cex = .9)


legend(-.2, -3.35, c("Without Identity Politics Loss Narrative", "With Identity Politics Loss Narrative ", expression(paste("Significant Effect, ", italic("p"), " < .10 "))), col = c("black"), pch = c(21, 16, 8), xpd = NA)


plotCI(blackmen.predprob$point, blackmen.predprob$yvalue, xlab = "", ui = blackmen.predprob$upper, li = blackmen.predprob$lower, sfrac = .002, err = "x", lwd = 1, gap = 0, pch = c(21, 16, 21, 16, 21, 8, 21, 16, 21, 16, 21, 16, 21, 8, 21, 8, 21, 8), yaxt = "n", ylab = "", main = "Black Men", xlim = c(-.31, .31), ylim = c(.5, 35.3), cex.axis = .8, cex.main = 1)
abline(v = 0, lty = 2)

lines(c(blackmen.predprob$upper[6] + .05, blackmen.predprob$upper[6] + .05), c(blackmen.predprob$yvalue[6], blackmen.predprob$yvalue[5]))
lines(c(blackmen.predprob$upper[6] + .05, blackmen.predprob$upper[6] + .02), c(blackmen.predprob$yvalue[6], blackmen.predprob$yvalue[6]))
lines(c(blackmen.predprob$upper[6] + .05, blackmen.predprob$upper[6] + .02), c(blackmen.predprob$yvalue[5], blackmen.predprob$yvalue[5]))
text(blackmen.predprob$upper[6] + .13, 24.5, "+.07", cex = .9)

lines(c(blackmen.predprob$upper[13] + .05, blackmen.predprob$upper[13] + .05), c(blackmen.predprob$yvalue[13], blackmen.predprob$yvalue[14]))
lines(c(blackmen.predprob$upper[13] + .05, blackmen.predprob$upper[13] + .02), c(blackmen.predprob$yvalue[13], blackmen.predprob$yvalue[13]))
lines(c(blackmen.predprob$upper[13] + .05, blackmen.predprob$upper[13] + .02), c(blackmen.predprob$yvalue[14], blackmen.predprob$yvalue[14]))
text(blackmen.predprob$upper[13] + .12, 7.7, "-.07", cex = .9)

lines(c(blackmen.predprob$upper[15] + .05, blackmen.predprob$upper[15] + .05), c(blackmen.predprob$yvalue[15], blackmen.predprob$yvalue[16]))
lines(c(blackmen.predprob$upper[15] + .05, blackmen.predprob$upper[15] + .02), c(blackmen.predprob$yvalue[15], blackmen.predprob$yvalue[15]))
lines(c(blackmen.predprob$upper[15] + .05, blackmen.predprob$upper[15] + .02), c(blackmen.predprob$yvalue[16], blackmen.predprob$yvalue[16]))
text(blackmen.predprob$upper[15] + .12, 4.5, "-.07", cex = .9)

lines(c(blackmen.predprob$upper[17] + .05, blackmen.predprob$upper[17] + .05), c(blackmen.predprob$yvalue[17], blackmen.predprob$yvalue[18]))
lines(c(blackmen.predprob$upper[17] + .05, blackmen.predprob$upper[17] + .02), c(blackmen.predprob$yvalue[17], blackmen.predprob$yvalue[17]))
lines(c(blackmen.predprob$upper[17] + .05, blackmen.predprob$upper[17] + .02), c(blackmen.predprob$yvalue[18], blackmen.predprob$yvalue[18]))
text(blackmen.predprob$upper[17] + .12, 1.3, "-.08", cex = .9)


plotCI(blackwomen.predprob$point, blackwomen.predprob$yvalue, xlab = "", ui = blackwomen.predprob$upper, li = blackwomen.predprob$lower, sfrac = .002, err = "x", lwd = 1, gap = 0, pch = c(21, 16), yaxt = "n", ylab = "", main = "Black Women", xlim = c(-.31, .31), ylim = c(.5, 35.3), cex.axis = .8, cex.main = 1)
abline(v = 0, lty = 2)


#####DEMOGRAPHIC CHARACTERISTICS (OA-C)#####

standard.ww <- subset(standard.w,female == 1)
standard.wm <- subset(standard.w,female == 0)
standard.bw <- subset(standard.b,female == 1)
standard.bm <- subset(standard.b,female == 0)

age.whitewomen <- mean(standard.ww$age)
print(age.whitewomen)

income.whitewomen <- median(standard.ww$income)
print(income.whitewomen)

education.whitewomen <- median(standard.ww$education)
print(education.whitewomen)

ideology.whitewomen <- mean(na.omit(standard.ww$ideology))
print(ideology.whitewomen)

married.whitewomen <- mean(standard.ww$married)
print(married.whitewomen)



age.whitemen <- mean(standard.wm$age)
print(age.whitemen)

income.whitemen <- median(standard.wm$income)
print(income.whitemen)

education.whitemen <- median(standard.wm$education)
print(education.whitemen)

ideology.whitemen <- mean(na.omit(standard.wm$ideology))
print(ideology.whitemen)

married.whitemen <- mean(standard.wm$married)
print(married.whitemen)



age.blackwomen <- mean(standard.bw$age)
print(age.blackwomen)

income.blackwomen <- median(na.omit(standard.bw$income))
print(income.blackwomen)

education.blackwomen <- median(standard.bw$education)
print(education.blackwomen)

ideology.blackwomen <- mean(na.omit(standard.bw$ideology))
print(ideology.blackwomen)

married.blackwomen <- mean(standard.bw$married)
print(married.blackwomen)



age.blackmen <- mean(standard.bm$age)
print(age.blackmen)

income.blackmen <- median(standard.bm$income)
print(income.blackmen)

education.blackmen <- median(na.omit(standard.bm$education))
print(education.blackmen)

ideology.blackmen <- mean(na.omit(standard.bm$ideology))
print(ideology.blackmen)

married.blackmen <- mean(standard.bm$married)
print(married.blackmen)


#####BALANCE TESTING (OA-D)#####

m.out1 <- MatchBalance(treatment ~ age + income + education + married + ideology + female, data = standard.w, nboots = 100000, na.omit) #May differ slightly due to simulation

m.out2 <- MatchBalance(treatment ~ age + income + education + married + ideology + female, data = standard.b, nboots = 100000, na.omit)





#####PROBABILITY BY CANDIDATE PROFILE (OA-G & OA-H)#####

#m1 is all
#m2 is men
#m3 is women

#Generate the predicted probabilities, view them in excel, manually input results into figure

#Get predicted probabilities for all candidates

allcombos.control <- expand.grid(treatment = 0, candgender = c(0, 1), candrace = c(1, 2, 3, 4), candideol = c(1, 2, 3), candpolicy = c(1, 2, 3, 4))
head(allcombos.control)
nrow(allcombos.control)

allcombos.treatment <- expand.grid(treatment = 1, candgender = c(0, 1), candrace = c(1, 2, 3, 4), candideol = c(1, 2, 3), candpolicy = c(1, 2, 3, 4))
head(allcombos.treatment)
nrow(allcombos.treatment)


#Probabilities for white men

control.all <- data.frame(matrix(ncol = 7, nrow = 0))
x <- c("treatment", "candgender", "candrace", "candideol", "candpolicy", "prob", "se")
colnames(control.all) <- x

treatment.all <- data.frame(matrix(ncol = 7, nrow = 0))
x <- c("treatment", "candgender", "candrace", "candideol", "candpolicy", "prob", "se")
colnames(treatment.all) <- x


m1 <- zelig(candvote ~ candgender*treatment + candrace*treatment + candideol*treatment + candpolicy*treatment, model = "logit.gee", id = "data.id", data = men.white, cite = FALSE)
summary(m1)


i <- 1
while (i <= 96){
  newdata <- setx(m1, treatment = allcombos.control[i, 1], candgender = allcombos.control[i, 2], candrace = allcombos.control[i, 3], candideol = allcombos.control[i, 4], candpolicy = allcombos.control[i, 5])
  predprob1 <- sim(m1, x = newdata, num = 100000)
  predprob2 <- zelig_qi_to_df((predprob1))
  prob <- mean(predprob2$expected_value)
  se <- sd(predprob2$expected_value)
  predprob.all <- cbind(allcombos.control[i, ], prob, se)
  control.all <- rbind(control.all, predprob.all)
  i <- i + 1
}


control.whitemen <- control.all
write.csv(control.whitemen, "White Men Control Probabilities.csv")



i <- 1
while (i <= 96){
  newdata <- setx(m1, treatment = allcombos.treatment[i, 1], candgender = allcombos.treatment[i, 2], candrace = allcombos.treatment[i, 3], candideol = allcombos.treatment[i, 4], candpolicy = allcombos.treatment[i, 5])
  predprob1 <- sim(m1, x = newdata, num = 100000)
  predprob2 <- zelig_qi_to_df((predprob1))
  prob <- mean(predprob2$expected_value)
  se <- sd(predprob2$expected_value)
  predprob.all <- cbind(allcombos.treatment[i, ], prob, se)
  treatment.all <- rbind(treatment.all, predprob.all)
  i <- i + 1
}


treatment.whitemen <- treatment.all
write.csv(treatment.whitemen, "White Men Treatment Probabilities.csv")





#Probabilities for white women

control.all <- data.frame(matrix(ncol = 7, nrow = 0))
x <- c("treatment", "candgender", "candrace", "candideol", "candpolicy", "prob", "se")
colnames(control.all) <- x

treatment.all <- data.frame(matrix(ncol = 7, nrow = 0))
x <- c("treatment", "candgender", "candrace", "candideol", "candpolicy", "prob", "se")
colnames(treatment.all) <- x


m1 <- zelig(candvote ~ candgender*treatment + candrace*treatment + candideol*treatment + candpolicy*treatment, model = "logit.gee", id = "data.id", data = women.white, cite = FALSE)
summary(m1)


i <- 1
while (i <= 96){
  newdata <- setx(m1, treatment = allcombos.control[i, 1], candgender = allcombos.control[i, 2], candrace = allcombos.control[i, 3], candideol = allcombos.control[i, 4], candpolicy = allcombos.control[i, 5])
  predprob1 <- sim(m1, x = newdata, num = 100000)
  predprob2 <- zelig_qi_to_df((predprob1))
  prob <- mean(predprob2$expected_value)
  se <- sd(predprob2$expected_value)
  predprob.all <- cbind(allcombos.control[i, ], prob, se)
  control.all <- rbind(control.all, predprob.all)
  i <- i + 1
}


control.whitewomen <- control.all
write.csv(control.whitewomen, "White Women Control Probabilities.csv")



i <- 1
while (i <= 96){
  newdata <- setx(m1, treatment = allcombos.treatment[i, 1], candgender = allcombos.treatment[i, 2], candrace = allcombos.treatment[i, 3], candideol = allcombos.treatment[i, 4], candpolicy = allcombos.treatment[i, 5])
  predprob1 <- sim(m1, x = newdata, num = 100000)
  predprob2 <- zelig_qi_to_df((predprob1))
  prob <- mean(predprob2$expected_value)
  se <- sd(predprob2$expected_value)
  predprob.all <- cbind(allcombos.treatment[i, ], prob, se)
  treatment.all <- rbind(treatment.all, predprob.all)
  i <- i + 1
}


treatment.whitewomen <- treatment.all
write.csv(treatment.whitewomen, "White Women Treatment Probabilities.csv")





#Probabilities for black men

control.all <- data.frame(matrix(ncol = 7, nrow = 0))
x <- c("treatment", "candgender", "candrace", "candideol", "candpolicy", "prob", "se")
colnames(control.all) <- x

treatment.all <- data.frame(matrix(ncol = 7, nrow = 0))
x <- c("treatment", "candgender", "candrace", "candideol", "candpolicy", "prob", "se")
colnames(treatment.all) <- x


m1 <- zelig(candvote ~ candgender*treatment + candrace*treatment + candideol*treatment + candpolicy*treatment, model = "logit.gee", id = "data.id", data = men.black, cite = FALSE)
summary(m1)


i <- 1
while (i <= 96){
  newdata <- setx(m1, treatment = allcombos.control[i, 1], candgender = allcombos.control[i, 2], candrace = allcombos.control[i, 3], candideol = allcombos.control[i, 4], candpolicy = allcombos.control[i, 5])
  predprob1 <- sim(m1, x = newdata, num = 100000)
  predprob2 <- zelig_qi_to_df((predprob1))
  prob <- mean(predprob2$expected_value)
  se <- sd(predprob2$expected_value)
  predprob.all <- cbind(allcombos.control[i, ], prob, se)
  control.all <- rbind(control.all, predprob.all)
  i <- i + 1
}


control.blackmen <- control.all
write.csv(control.blackmen, "Black Men Control Probabilities.csv")



i <- 1
while (i <= 96){
  newdata <- setx(m1, treatment = allcombos.treatment[i, 1], candgender = allcombos.treatment[i, 2], candrace = allcombos.treatment[i, 3], candideol = allcombos.treatment[i, 4], candpolicy = allcombos.treatment[i, 5])
  predprob1 <- sim(m1, x = newdata, num = 100000)
  predprob2 <- zelig_qi_to_df((predprob1))
  prob <- mean(predprob2$expected_value)
  se <- sd(predprob2$expected_value)
  predprob.all <- cbind(allcombos.treatment[i, ], prob, se)
  treatment.all <- rbind(treatment.all, predprob.all)
  i <- i + 1
}


treatment.blackmen<- treatment.all
write.csv(treatment.blackmen, "Black Men Treatment Probabilities.csv")





#Probabilities for black women

control.all <- data.frame(matrix(ncol = 7, nrow = 0))
x <- c("treatment", "candgender", "candrace", "candideol", "candpolicy", "prob", "se")
colnames(control.all) <- x

treatment.all <- data.frame(matrix(ncol = 7, nrow = 0))
x <- c("treatment", "candgender", "candrace", "candideol", "candpolicy", "prob", "se")
colnames(treatment.all) <- x


m1 <- zelig(candvote ~ candgender*treatment + candrace*treatment + candideol*treatment + candpolicy*treatment, model = "logit.gee", id = "data.id", data = women.black, cite = FALSE)
summary(m1)


i <- 1
while (i <= 96){
  newdata <- setx(m1, treatment = allcombos.control[i, 1], candgender = allcombos.control[i, 2], candrace = allcombos.control[i, 3], candideol = allcombos.control[i, 4], candpolicy = allcombos.control[i, 5])
  predprob1 <- sim(m1, x = newdata, num = 100000)
  predprob2 <- zelig_qi_to_df((predprob1))
  prob <- mean(predprob2$expected_value)
  se <- sd(predprob2$expected_value)
  predprob.all <- cbind(allcombos.control[i, ], prob, se)
  control.all <- rbind(control.all, predprob.all)
  i <- i + 1
}


control.blackwomen <- control.all
write.csv(control.blackwomen, "Black Women Control Probabilities.csv")



i <- 1
while (i <= 96){
  newdata <- setx(m1, treatment = allcombos.treatment[i, 1], candgender = allcombos.treatment[i, 2], candrace = allcombos.treatment[i, 3], candideol = allcombos.treatment[i, 4], candpolicy = allcombos.treatment[i, 5])
  predprob1 <- sim(m1, x = newdata, num = 100000)
  predprob2 <- zelig_qi_to_df((predprob1))
  prob <- mean(predprob2$expected_value)
  se <- sd(predprob2$expected_value)
  predprob.all <- cbind(allcombos.treatment[i, ], prob, se)
  treatment.all <- rbind(treatment.all, predprob.all)
  i <- i + 1
}


treatment.blackwomen <- treatment.all
write.csv(treatment.blackwomen, "Black Women Treatment Probabilities.csv")



#Probability of Winning by Candidate Traits

#Derived from analysis in Excel; results may vary slightly due variation in simulation

simulation.probs <- read.csv("SimulationProbs.csv")
head(simulation.probs)

prob.whitemen.control <- simulation.probs$white.men.control
prob.whitemen.treatment <- c(simulation.probs$white.men.treatment)
prob.whitemen <- rbind(prob.whitemen.control, prob.whitemen.treatment)
head(prob.whitemen)


prob.whitewomen.control <- simulation.probs$white.women.control
prob.whitewomen.treatment <- c(simulation.probs$white.women.treatment)
prob.whitewomen <- rbind(prob.whitewomen.control, prob.whitewomen.treatment)
head(prob.whitewomen)


prob.blackmen.control <- simulation.probs$black.men.control
prob.blackmen.treatment <- c(simulation.probs$black.men.treatment)
prob.blackmen <- rbind(prob.blackmen.control, prob.blackmen.treatment)
head(prob.blackmen)


prob.blackwomen.control <- simulation.probs$black.women.control
prob.blackwomen.treatment <- c(simulation.probs$black.women.treatment)
prob.blackwomen <- rbind(prob.blackwomen.control, prob.blackwomen.treatment)
head(prob.blackwomen)


#Dimensions are 6.5x6.5
par(family = "serif")
par(omi = c(.75, 0, 0, 0))
par(mar = c(2, 3.5, 2.5, .25))
par(mfrow = c(4, 1))
groups <- c("Man", "Woman", "White", "Asian", "Latino", "Black", "Moderately\nLiberal", "Slightly\nLiberal", "Very\nLiberal", "Good\nJobs", "Justice\nReform", "Fair\nEmployment", "Fair\nSentencing")

x <- barplot(prob.whitemen, ylim = c(0, 1.1), xaxt = "n", col = c("grey80", "grey40"), beside = TRUE, space =c(0, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0), cex.axis = 1)
box()
xspot <- x[1,] + (x[2,] - x[1,])/2
text(17.5, 1.18, xpd = TRUE, "White Democratic Men", cex = 1.2)
text(xspot, -.11, xpd = T, c("Man\n ", "Woman\n ", "White\n ", "Asian\n ", "Latino\n ", "Black\n ", "Moderate\nLiberal", "Slight\nLiberal", "Very\nLiberal", "Good\nJobs", "Justice\nReform", "Fair\nEmploy", "Fair\nSentence"), cex = 1)
text(x[ ,], prob.whitemen + .05, rd(prob.whitemen, digits = 2), cex = .9)


x <- barplot(prob.whitewomen, ylim = c(0, 1.1), xaxt = "n", col = c("grey80", "grey40"), beside = TRUE, space =c(0, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0), cex.axis = 1)
box()
xspot <- x[1,] + (x[2,] - x[1,])/2
text(17.5, 1.18, xpd = TRUE, "White Democratic Women", cex = 1.2)
text(xspot, -.11, xpd = T, c("Man\n ", "Woman\n ", "White\n ", "Asian\n ", "Latino\n ", "Black\n ", "Moderate\nLiberal", "Slight\nLiberal", "Very\nLiberal", "Good\nJobs", "Justice\nReform", "Fair\nEmploy", "Fair\nSentence"), cex = 1)
text(x[ ,], prob.whitewomen + .05, rd(prob.whitewomen, digits = 2), cex = .9)


x <- barplot(prob.blackmen, ylim = c(0, 1.1), xaxt = "n", col = c("grey80", "grey40"), beside = TRUE, space =c(0, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0), cex.axis = 1)
box()
xspot <- x[1,] + (x[2,] - x[1,])/2
text(17.5, 1.18, xpd = TRUE, "Black Democratic Men", cex = 1.2)
text(xspot, -.11, xpd = T, c("Man\n ", "Woman\n ", "White\n ", "Asian\n ", "Latino\n ", "Black\n ", "Moderate\nLiberal", "Slight\nLiberal", "Very\nLiberal", "Good\nJobs", "Justice\nReform", "Fair\nEmploy", "Fair\nSentence"), cex = 1)
text(x[ ,], prob.blackmen + .05, rd(prob.blackmen, digits = 2), cex = .9)


x <- barplot(prob.blackwomen, ylim = c(0, 1.1), xaxt = "n", col = c("grey80", "grey40"), beside = TRUE, space =c(0, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0, .75, 0), cex.axis = 1)
box()
xspot <- x[1,] + (x[2,] - x[1,])/2
text(17.5, 1.18, xpd = TRUE, "Black Democratic Women", cex = 1.2)
text(xspot, -.11, xpd = T, c("Man\n ", "Woman\n ", "White\n ", "Asian\n ", "Latino\n ", "Black\n ", "Moderate\nLiberal", "Slight\nLiberal", "Very\nLiberal", "Good\nJobs", "Justice\nReform", "Fair\nEmploy", "Fair\nSentence"), cex = 1)
text(x[ ,], prob.blackwomen + .05, rd(prob.blackwomen, digits = 2), cex = .9)

par(fig=c(0, 1, 0, 1), oma=c(0, 0, 0, 0), mar=c(0, 0, 0, 0), new=TRUE)
plot(0, 0, type='n', bty='n', xaxt='n', yaxt='n')
legend(-.926, -.95, border = "black", fill = c("grey80", "grey40"), c("Without Loss Narrative", "With Loss Narrative"))








#####REGRESSIONS FOR CANDIDATE VOTE CHOICE BY RACE AND GENDER (OA-I)#####

#Preliminary Regressions
wmc <- zelig(candvote ~ candgender + candrace + candideol + candpolicy, model = "logit.gee", id = "data.id", data = men.white.control, cite = FALSE) #White men in control
summary(wmc) #Ignore any error messages

wmt <- zelig(candvote ~ candgender + candrace + candideol + candpolicy, model = "logit.gee", id = "data.id", data = men.white.treat, cite = FALSE) #White men in treatment
summary(wmt)

wwc <- zelig(candvote ~ candgender + candrace + candideol + candpolicy, model = "logit.gee", id = "data.id", data = women.white.control, cite = FALSE) #White women in control
summary(wwc)

wwt <- zelig(candvote ~ candgender + candrace + candideol + candpolicy, model = "logit.gee", id = "data.id", data = women.white.treat, cite = FALSE) #White women in treatment
summary(wwt)


bmc <- zelig(candvote ~ candgender + candrace + candideol + candpolicy, model = "logit.gee", id = "data.id", data = men.black.control, cite = FALSE) #Black men in control
summary(bmc)

bmt <- zelig(candvote ~ candgender + candrace + candideol + candpolicy, model = "logit.gee", id = "data.id", data = men.black.treat, cite = FALSE) #Black men in treatment
summary(bmt)

bwc <- zelig(candvote ~ candgender + candrace + candideol + candpolicy, model = "logit.gee", id = "data.id", data = women.black.control, cite = FALSE) #Black women in control
summary(bwc)

bwt <- zelig(candvote ~ candgender + candrace + candideol + candpolicy, model = "logit.gee", id = "data.id", data = women.black.treat, cite = FALSE) #Black women in treatment
summary(bwt)


#Main Regressions
logit.white.men <- zelig(candvote ~ candgender*treatment + candrace*treatment + candideol*treatment + candpolicy*treatment, model = "logit.gee", id = "data.id", data = men.white, cite = FALSE) #White men
summary(logit.white.men)
nrow(na.omit(men.white))


logit.white.women <- zelig(candvote ~ candgender*treatment + candrace*treatment + candideol*treatment + candpolicy*treatment, model = "logit.gee", id = "data.id", data = women.white, cite = FALSE) #White women
summary(logit.white.women)
nrow(na.omit(women.white))


logit.black.men <- zelig(candvote ~ candgender*treatment + candrace*treatment + candideol*treatment + candpolicy*treatment, model = "logit.gee", id = "data.id", data = men.black, cite = FALSE)
summary(logit.black.men)
nrow(na.omit(men.black))


logit.black.women <- zelig(candvote ~ candgender*treatment + candrace*treatment + candideol*treatment + candpolicy*treatment, model = "logit.gee", id = "data.id", data = women.black, cite = FALSE)
summary(logit.black.women)
nrow(na.omit(women.black))



#####REGRESSIONS FOR CANDIDATE VOTE CHOICE, MAIN EFFECTS CONTROL (OA-J)#####

#Organize Data

smallstackedb <- subset(stacked.b, select = c(id, candvote, treatment, candgender, candrace, candideology, candpolicy, female), treatment == 0) #Create a subset that includes only the variables needed for the regressions
head(smallstackedb)
nrow(smallstackedb)

smallstackedb$candrace <- as.factor(smallstackedb$candrace) #Create factor variables, white, moderate, good jobs, male as the baseline
smallstackedb$candpolicy <- as.factor(smallstackedb$candpolicy)
smallstackedb$candideology <- as.factor(smallstackedb$candideology)

women.black <- subset(smallstackedb, female == 1)
men.black <- subset(smallstackedb, female == 0)

smallstackedw <- subset(stacked.w, select = c(id, candvote, treatment, candgender, candrace, candideology, candpolicy, female), treatment == 0) #Create a subset that includes only the variables needed for the regressions
head(smallstackedw)
nrow(smallstackedw)

smallstackedw$candrace <- as.factor(smallstackedw$candrace) #Create factor variables, white, moderate, good jobs, male as the baseline
smallstackedw$candpolicy <- as.factor(smallstackedw$candpolicy)
smallstackedw$candideology <- as.factor(smallstackedw$candideology)

women.white <- subset(smallstackedw, female == 1)
men.white <- subset(smallstackedw, female == 0)

women.black.control <- subset(women.black, treatment == 0)
women.black.treat <- subset(women.black, treatment == 1)
men.black.control <- subset(men.black, treatment == 0)
men.black.treat <- subset(men.black, treatment == 1)

women.white.control <- subset(women.white, treatment == 0)
women.white.treat <- subset(women.white, treatment == 1)
men.white.control <- subset(men.white, treatment == 0)
men.white.treat <- subset(men.white, treatment == 1)

mwc <- zelig(candvote ~ candgender + candrace + candideology + candpolicy, model = "logit.gee", id = "id", data = men.white.control, cite = FALSE)
summary(mwc)
nrow(na.omit(men.white.control))

wwc <- zelig(candvote ~ candgender + candrace + candideology + candpolicy, model = "logit.gee", id = "id", data = women.white.control, cite = FALSE)
summary(wwc)
nrow(na.omit(women.white.control))

mbc <- zelig(candvote ~ candgender + candrace + candideology + candpolicy, model = "logit.gee", id = "id", data = men.black.control, cite = FALSE)
summary(mbc)
nrow(na.omit(men.black.control))

wbc <- zelig(candvote ~ candgender + candrace + candideology + candpolicy, model = "logit.gee", id = "id", data = women.black.control, cite = FALSE)
summary(wbc)
nrow(na.omit(women.black.control))



#####REGRESSIONS FOR FIRST HEAT ROBUSTNESS CHECK (OA-K)#####

smallstacked.robust <- subset(stacked, heat == 1, select = c(responsenum, heat, candvote, treatment, candgender, candrace, candideology, candideolorder, candpolicy, female)) #Create a subset that includes only the variables needed for the regressions
head(smallstacked.robust)
nrow(smallstacked.robust)

smallstacked.robust$candrace <- as.factor(smallstacked.robust$candrace) #Create factor variables, white, moderate, good jobs, male as the baseline
smallstacked.robust$candpolicy <- as.factor(smallstacked.robust$candpolicy)
smallstacked.robust$candideolorder <- as.factor(smallstacked.robust$candideolorder)

women.small <- subset(smallstacked.robust, female == 1) #Create subsets by respondent gender and condition
men.small <- subset(smallstacked.robust, female == 0)
all.small <- smallstacked.robust
women.smalltreat <- subset(smallstacked.robust, treatment == 1)
women.smallcont <- subset(smallstacked.robust, treatment == 0)
men.smalltreat <- subset(smallstacked.robust, treatment == 1)
men.smallcont <- subset(smallstacked.robust, treatment == 0)
all.smallcont <- subset(smallstacked.robust, treatment == 0)
all.smalltreat <- subset(smallstacked.robust, treatment == 1)


m1 <- zelig(candvote ~ candgender*treatment + candrace*treatment + candideolorder*treatment + candpolicy*treatment, model = "logit.gee", id = "responsenum", data = all.small, cite = FALSE)
summary(m1)
nrow(na.omit(all.small))

m2 <- zelig(candvote ~ candgender*treatment + candrace*treatment + candideolorder*treatment + candpolicy*treatment, model = "logit.gee", id = "responsenum", data = men.small, cite = FALSE)
summary(m2)
nrow(na.omit(men.small))

m3 <- zelig(candvote ~ candgender*treatment + candrace*treatment + candideolorder*treatment + candpolicy*treatment, model = "logit.gee", id = "responsenum", data = women.small, cite = FALSE)
summary(m3)
nrow(na.omit(women.small))

#####MARGINAL PROBABILITY OF VOTE CHOICE IN FIRST HEAT (OA-K)#####

#m1 is all
#m2 is men
#m3 is women

#Get predicted probabilities for all respondents

#All respondents, gender
predprob1 <- setx(m1, treatment = 0, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1) #Get predicted probabilities for overall gender
predprob2 <- setx(m1, treatment = 0, candgender = 1, candrace = 1, candideolorder = 0, candpolicy = 1)
predprob3 <- setx(m1, treatment = 1, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1)
predprob4 <- setx(m1, treatment = 1, candgender = 1, candrace = 1, candideolorder = 0, candpolicy = 1)

gender.controlfull <- sim(m1, x = predprob1, x1 = predprob2, num = 100000)
summary(gender.controlfull)

gender.treatfull <- sim(m1, x = predprob3, x1 = predprob4, num = 1000)
summary(gender.treatfull)

genderdiff.full <- c(.08, .072)
genderdiff.fullbar <- genderdiff.full - c(.021, .017)


#All respondents, race
predprob1 <- setx(m1, treatment = 0, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1) #Get predicted probabilities for overall race
predprob2 <- setx(m1, treatment = 0, candgender = 0, candrace = 2, candideolorder = 0, candpolicy = 1)
predprob3 <- setx(m1, treatment = 0, candgender = 0, candrace = 3, candideolorder = 0, candpolicy = 1)
predprob4 <- setx(m1, treatment = 0, candgender = 0, candrace = 4, candideolorder = 0, candpolicy = 1)
predprob5 <- setx(m1, treatment = 1, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1)
predprob6 <- setx(m1, treatment = 1, candgender = 0, candrace = 2, candideolorder = 0, candpolicy = 1)
predprob7 <- setx(m1, treatment = 1, candgender = 0, candrace = 3, candideolorder = 0, candpolicy = 1)
predprob8 <- setx(m1, treatment = 1, candgender = 0, candrace = 4, candideolorder = 0, candpolicy = 1)

raceasian.controlfull <- sim(m1, x = predprob1, x1 = predprob2, num = 100000)
summary(raceasian.controlfull)

raceasian.treatfull <- sim(m1, x = predprob5, x1 = predprob6, num = 100000)
summary(raceasian.treatfull)

racehispanic.controlfull <- sim(m1, x = predprob1, x1 = predprob3, num = 100000)
summary(racehispanic.controlfull)

racehispanic.treatfull <- sim(m1, x = predprob5, x1 = predprob7, num = 100000)
summary(racehispanic.treatfull)

raceblack.controlfull <- sim(m1, x = predprob1, x1 = predprob4, num = 100000)
summary(raceblack.controlfull)

raceblack.treatfull <- sim(m1, x = predprob5, x1 = predprob8, num = 100000)
summary(raceblack.treatfull)


racediff.full <- c(-0.08, -.062, -.133, -.048, -.028, .003)
racediff.fullbar <- racediff.full- c(-.174, -.147, -.222, -.131, -.113, -.079)



#All respondents, ideology
predprob1 <- setx(m1, treatment = 0, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1) #Get predicted probabilities for overall ideology
predprob2 <- setx(m1, treatment = 0, candgender = 0, candrace = 1, candideolorder = 1, candpolicy = 1)
predprob3 <- setx(m1, treatment = 0, candgender = 0, candrace = 1, candideolorder = 2, candpolicy = 1)
predprob4 <- setx(m1, treatment = 1, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1)
predprob5 <- setx(m1, treatment = 1, candgender = 0, candrace = 1, candideolorder = 1, candpolicy = 1)
predprob6 <- setx(m1, treatment = 1, candgender = 0, candrace = 1, candideolorder = 2, candpolicy = 1)

ideologycon.controlfull <- sim(m1, x = predprob1, x1 = predprob2, num = 100000)
summary(ideologycon.controlfull)

ideologycon.treatfull <- sim(m1, x = predprob4, x1 = predprob5, num = 100000)
summary(ideologycon.treatfull)

ideologylib.controlfull <- sim(m1, x = predprob1, x1 = predprob3, num = 100000)
summary(ideologylib.controlfull)

ideologylib.treatfull <- sim(m1, x = predprob4, x1 = predprob6, num = 100000)
summary(ideologylib.treatfull)


ideologydiff.full <- c(-.082, -.161, -.113, -.128)
ideologydiff.fullbar <- ideologydiff.full - c(-.163, -.239, -.192, -.206)



#All respondents, policy
predprob1 <- setx(m1, treatment = 0, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1)
predprob2 <- setx(m1, treatment = 0, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 3)
predprob3 <- setx(m1, treatment = 0, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 2)
predprob4 <- setx(m1, treatment = 0, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 4)
predprob5 <- setx(m1, treatment = 1, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1)
predprob6 <- setx(m1, treatment = 1, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 3)
predprob7 <- setx(m1, treatment = 1, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 2)
predprob8 <- setx(m1, treatment = 1, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 4)


policyreform.controlfull <- sim(m1, x = predprob1, x1 = predprob2, num = 100000)
summary(policyreform.controlfull)

policyreform.treatfull <- sim(m1, x = predprob5, x1 = predprob6, num = 100000)
summary(policyreform.treatfull)

policyfairjob.controlfull <- sim(m1, x = predprob1, x1 = predprob3, num = 100000)
summary(policyfairjob.controlfull)

policyfairjob.treatfull <- sim(m1, x = predprob5, x1 = predprob7, num = 100000)
summary(policyfairjob.treatfull)

policyfairprison.controlfull <- sim(m1, x = predprob1, x1 = predprob4, num = 100000)
summary(policyfairprison.controlfull)

policyfairprison.treatfull <- sim(m1, x = predprob5, x1 = predprob8, num = 100000)
summary(policyfairprison.treatfull)


policydiff.full <- c(-.102, -.203, -.107, -.127, -.112, -.208)
policydiff.fullbar <- policydiff.full - c(-.192, -.298, -.203, -.217, -.203, -.301)


all.predprob <- c(policydiff.full, ideologydiff.full, racediff.full, genderdiff.full)
all.ci <- c(policydiff.fullbar, ideologydiff.fullbar, racediff.fullbar, genderdiff.fullbar)





#Get predicted probabilities for men respondents

#men respondents, gender
predprob1 <- setx(m2, treatment = 0, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1) #Get predicted probabilities for overmen gender
predprob2 <- setx(m2, treatment = 0, candgender = 1, candrace = 1, candideolorder = 0, candpolicy = 1)
predprob3 <- setx(m2, treatment = 1, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1)
predprob4 <- setx(m2, treatment = 1, candgender = 1, candrace = 1, candideolorder = 0, candpolicy = 1)

gender.controlmen <- sim(m2, x = predprob1, x1 = predprob2, num = 100000)
summary(gender.controlmen)

gender.treatmen <- sim(m2, x = predprob3, x1 = predprob4, num = 1000)
summary(gender.treatmen)

genderdiff.men <- c(.010, .061)
genderdiff.menbar <- genderdiff.men - c(-.056, -.023)


#men respondents, race
predprob1 <- setx(m2, treatment = 0, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1) #Get predicted probabilities for overmen race
predprob2 <- setx(m2, treatment = 0, candgender = 0, candrace = 2, candideolorder = 0, candpolicy = 1)
predprob3 <- setx(m2, treatment = 0, candgender = 0, candrace = 3, candideolorder = 0, candpolicy = 1)
predprob4 <- setx(m2, treatment = 0, candgender = 0, candrace = 4, candideolorder = 0, candpolicy = 1)
predprob5 <- setx(m2, treatment = 1, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1)
predprob6 <- setx(m2, treatment = 1, candgender = 0, candrace = 2, candideolorder = 0, candpolicy = 1)
predprob7 <- setx(m2, treatment = 1, candgender = 0, candrace = 3, candideolorder = 0, candpolicy = 1)
predprob8 <- setx(m2, treatment = 1, candgender = 0, candrace = 4, candideolorder = 0, candpolicy = 1)

raceasian.controlmen <- sim(m2, x = predprob1, x1 = predprob2, num = 100000)
summary(raceasian.controlmen)

raceasian.treatmen <- sim(m2, x = predprob5, x1 = predprob6, num = 100000)
summary(raceasian.treatmen)

racehispanic.controlmen <- sim(m2, x = predprob1, x1 = predprob3, num = 100000)
summary(racehispanic.controlmen)

racehispanic.treatmen <- sim(m2, x = predprob5, x1 = predprob7, num = 100000)
summary(racehispanic.treatmen)

raceblack.controlmen <- sim(m2, x = predprob1, x1 = predprob4, num = 100000)
summary(raceblack.controlmen)

raceblack.treatmen <- sim(m2, x = predprob5, x1 = predprob8, num = 100000)
summary(raceblack.treatmen)


racediff.men <- c(-.115, -.102, -.181, -.075, -.116, -.041)
racediff.menbar <- racediff.men- c(-.231, -.221, -.301, -.201, -.221, -.160)



#men respondents, ideology
predprob1 <- setx(m2, treatment = 0, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1) #Get predicted probabilities for overmen ideology
predprob2 <- setx(m2, treatment = 0, candgender = 0, candrace = 1, candideolorder = 1, candpolicy = 1)
predprob3 <- setx(m2, treatment = 0, candgender = 0, candrace = 1, candideolorder = 2, candpolicy = 1)
predprob4 <- setx(m2, treatment = 1, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1)
predprob5 <- setx(m2, treatment = 1, candgender = 0, candrace = 1, candideolorder = 1, candpolicy = 1)
predprob6 <- setx(m2, treatment = 1, candgender = 0, candrace = 1, candideolorder = 2, candpolicy = 1)

ideologycon.controlmen <- sim(m2, x = predprob1, x1 = predprob2, num = 100000)
summary(ideologycon.controlmen)

ideologycon.treatmen <- sim(m2, x = predprob4, x1 = predprob5, num = 100000)
summary(ideologycon.treatmen)

ideologylib.controlmen <- sim(m2, x = predprob1, x1 = predprob3, num = 100000)
summary(ideologylib.controlmen)

ideologylib.treatmen <- sim(m2, x = predprob4, x1 = predprob6, num = 100000)
summary(ideologylib.treatmen)


ideologydiff.men <- c(-.055, -.113, -.093, -.067)
ideologydiff.menbar <- ideologydiff.men - c(-.143, -.221, -.186, -.175)



#men respondents, policy
predprob1 <- setx(m2, treatment = 0, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1)
predprob2 <- setx(m2, treatment = 0, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 3)
predprob3 <- setx(m2, treatment = 0, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 2)
predprob4 <- setx(m2, treatment = 0, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 4)
predprob5 <- setx(m2, treatment = 1, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1)
predprob6 <- setx(m2, treatment = 1, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 3)
predprob7 <- setx(m2, treatment = 1, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 2)
predprob8 <- setx(m2, treatment = 1, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 4)


policyreform.controlmen <- sim(m2, x = predprob1, x1 = predprob2, num = 100000)
summary(policyreform.controlmen)

policyreform.treatmen <- sim(m2, x = predprob5, x1 = predprob6, num = 100000)
summary(policyreform.treatmen)

policyfairjob.controlmen <- sim(m2, x = predprob1, x1 = predprob3, num = 100000)
summary(policyfairjob.controlmen)

policyfairjob.treatmen <- sim(m2, x = predprob5, x1 = predprob7, num = 100000)
summary(policyfairjob.treatmen)

policyfairprison.controlmen <- sim(m2, x = predprob1, x1 = predprob4, num = 100000)
summary(policyfairprison.controlmen)

policyfairprison.treatmen <- sim(m2, x = predprob5, x1 = predprob8, num = 100000)
summary(policyfairprison.treatmen)


policydiff.men <- c(-.175, -.119, -.154, -.094, -.144, -.188)
policydiff.menbar <- policydiff.men - c(-.281, -.246, -.286, -.224, -.259, -.317)


men.predprob <- c(policydiff.men, ideologydiff.men, racediff.men, genderdiff.men)
men.ci <- c(policydiff.menbar, ideologydiff.menbar, racediff.menbar, genderdiff.menbar)









#Get predicted probabilities for women respondents

#women respondents, gender
predprob1 <- setx(m3, treatment = 0, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1) #Get predicted probabilities for overwomen gender
predprob2 <- setx(m3, treatment = 0, candgender = 1, candrace = 1, candideolorder = 0, candpolicy = 1)
predprob3 <- setx(m3, treatment = 1, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1)
predprob4 <- setx(m3, treatment = 1, candgender = 1, candrace = 1, candideolorder = 0, candpolicy = 1)

gender.controlwomen <- sim(m3, x = predprob1, x1 = predprob2, num = 100000)
summary(gender.controlwomen)

gender.treatwomen <- sim(m3, x = predprob3, x1 = predprob4, num = 1000)
summary(gender.treatwomen)

genderdiff.women <- c(.17, .081)
genderdiff.womenbar <- genderdiff.women - c(.073, .012)


#women respondents, race
predprob1 <- setx(m3, treatment = 0, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1) #Get predicted probabilities for overwomen race
predprob2 <- setx(m3, treatment = 0, candgender = 0, candrace = 2, candideolorder = 0, candpolicy = 1)
predprob3 <- setx(m3, treatment = 0, candgender = 0, candrace = 3, candideolorder = 0, candpolicy = 1)
predprob4 <- setx(m3, treatment = 0, candgender = 0, candrace = 4, candideolorder = 0, candpolicy = 1)
predprob5 <- setx(m3, treatment = 1, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1)
predprob6 <- setx(m3, treatment = 1, candgender = 0, candrace = 2, candideolorder = 0, candpolicy = 1)
predprob7 <- setx(m3, treatment = 1, candgender = 0, candrace = 3, candideolorder = 0, candpolicy = 1)
predprob8 <- setx(m3, treatment = 1, candgender = 0, candrace = 4, candideolorder = 0, candpolicy = 1)

raceasian.controlwomen <- sim(m3, x = predprob1, x1 = predprob2, num = 100000)
summary(raceasian.controlwomen)

raceasian.treatwomen <- sim(m3, x = predprob5, x1 = predprob6, num = 100000)
summary(raceasian.treatwomen)

racehispanic.controlwomen <- sim(m3, x = predprob1, x1 = predprob3, num = 100000)
summary(racehispanic.controlwomen)

racehispanic.treatwomen <- sim(m3, x = predprob5, x1 = predprob7, num = 100000)
summary(racehispanic.treatwomen)

raceblack.controlwomen <- sim(m3, x = predprob1, x1 = predprob4, num = 100000)
summary(raceblack.controlwomen)

raceblack.treatwomen <- sim(m3, x = predprob5, x1 = predprob8, num = 100000)
summary(raceblack.treatwomen)


racediff.women <- c(-.013, -.026, -.073, -.033, .079, .047)
racediff.womenbar <- racediff.women- c(-.158, -.141, -.201, -.138, -.051, -.055)



#women respondents, ideology
predprob1 <- setx(m3, treatment = 0, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1) #Get predicted probabilities for overwomen ideology
predprob2 <- setx(m3, treatment = 0, candgender = 0, candrace = 1, candideolorder = 1, candpolicy = 1)
predprob3 <- setx(m3, treatment = 0, candgender = 0, candrace = 1, candideolorder = 2, candpolicy = 1)
predprob4 <- setx(m3, treatment = 1, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1)
predprob5 <- setx(m3, treatment = 1, candgender = 0, candrace = 1, candideolorder = 1, candpolicy = 1)
predprob6 <- setx(m3, treatment = 1, candgender = 0, candrace = 1, candideolorder = 2, candpolicy = 1)

ideologycon.controlwomen <- sim(m3, x = predprob1, x1 = predprob2, num = 100000)
summary(ideologycon.controlwomen)

ideologycon.treatwomen <- sim(m3, x = predprob4, x1 = predprob5, num = 100000)
summary(ideologycon.treatwomen)

ideologylib.controlwomen <- sim(m3, x = predprob1, x1 = predprob3, num = 100000)
summary(ideologylib.controlwomen)

ideologylib.treatwomen <- sim(m3, x = predprob4, x1 = predprob6, num = 100000)
summary(ideologylib.treatwomen)


ideologydiff.women <- c(-.091, -.211, -.122, -.186)
ideologydiff.womenbar <- ideologydiff.women - c(-.212, -.325, -.24, -.30)



#women respondents, policy
predprob1 <- setx(m3, treatment = 0, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1)
predprob2 <- setx(m3, treatment = 0, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 3)
predprob3 <- setx(m3, treatment = 0, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 2)
predprob4 <- setx(m3, treatment = 0, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 4)
predprob5 <- setx(m3, treatment = 1, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 1)
predprob6 <- setx(m3, treatment = 1, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 3)
predprob7 <- setx(m3, treatment = 1, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 2)
predprob8 <- setx(m3, treatment = 1, candgender = 0, candrace = 1, candideolorder = 0, candpolicy = 4)


policyreform.controlwomen <- sim(m3, x = predprob1, x1 = predprob2, num = 100000)
summary(policyreform.controlwomen)

policyreform.treatwomen <- sim(m3, x = predprob5, x1 = predprob6, num = 100000)
summary(policyreform.treatwomen)

policyfairjob.controlwomen <- sim(m3, x = predprob1, x1 = predprob3, num = 100000)
summary(policyfairjob.controlwomen)

policyfairjob.treatwomen <- sim(m3, x = predprob5, x1 = predprob7, num = 100000)
summary(policyfairjob.treatwomen)

policyfairprison.controlwomen <- sim(m3, x = predprob1, x1 = predprob4, num = 100000)
summary(policyfairprison.controlwomen)

policyfairprison.treatwomen <- sim(m3, x = predprob5, x1 = predprob8, num = 100000)
summary(policyfairprison.treatwomen)


policydiff.women <- c(-.004, -.307, -.037, -.16, -.043, -.225)
policydiff.womenbar <- policydiff.women - c(-.152, -.447, -.174, -.284, -.182, -.357)


women.predprob <- c(policydiff.women, ideologydiff.women, racediff.women, genderdiff.women)
women.ci <- c(policydiff.womenbar, ideologydiff.womenbar, racediff.womenbar, genderdiff.womenbar)


fullyvalue <- c(1, 1.5, 2.2, 4.2, 4.7, 5.4, 7.4, 7.9, 8.6,9.6, 12.6, 13.1, 13.8,15.8, 16.3, 17.0, 18.0,    21.0, 21.5, 22.2, 24.2, 24.7, 25.4,   27.4, 27.9, 28.6,   29.6,   32.6, 33.1, 33.8,   34.8)
pointyvalue <- c(1, 1.5, 4.2, 4.7, 7.4, 7.9, 12.6, 13.1, 15.8, 16.3, 21.0, 21.5, 24.2, 24.7, 27.4, 27.9, 32.6, 33.1)
groupyvalue <- c(2.2, 5.4, 8.6, 13.8, 17.0, 22.2, 25.4, 28.6, 33.8)
labyvalue <- c(9.6, 18.0, 29.6, 34.8)



par(mar=c(3, .45, 2, .45), oma = c(5, 12.0, 0, 0))
par(mfrow = c(1,3))
par(family = "serif")
plotCI(all.predprob, pointyvalue, xlab = "", uiw = all.ci, sfrac = .002, err = "x", lwd = 1, gap = 0, pch = c(21, 16), yaxt = "n", ylab = "", main = "All Respondents", xlim = c(-.50, .50), ylim = c(.5, 35.3), cex.axis = .8, cex.main = 1)
text(-.53, pointyvalue, pos = 2, c("Control", "Treatment", "Control", "Treatment","Control", "Treatment","Control", "Treatment","Control", "Treatment","Control", "Treatment") , srt = 360, xpd = NA)
text(-.53, groupyvalue, c("Fair Sentencing", "Fair Employment", "Justice Reform", "Very Liberal", "Slightly Liberal", "African American", "Hispanic", "Asian American", "Female"), pos = 2, font = 3, xpd = NA)
text(-.53, labyvalue, c("Policy (Good Jobs Baseline)", "Ideology (Moderate Baseline)", "Race (White Basline)", "Gender (Male Baseline)"), pos = 2, xpd = NA, font = 2)
abline(v = 0, lty = 2)



plotCI(men.predprob, pointyvalue, xlab = "", uiw = men.ci, sfrac = .002, err = "x", lwd = 1, gap = 0, pch = c(21, 16), yaxt = "n", ylab = "", main = "Men", xlim = c(-.50, .50), ylim = c(.5, 35.3), cex.axis = .8, cex.main = 1)
abline(v = 0, lty = 2)
legend(-.71, -3.35, c("With Identity Politics Loss Narrative", "Without Identity Politics Loss Narrative ", expression(paste("Significant Effect, ", italic("p"), " < .10 "))), col = c("black"), pch = c(16, 21, 8), xpd = NA)


plotCI(women.predprob, pointyvalue, xlab = "", uiw = women.ci, sfrac = .002, err = "x", lwd = 1, gap = 0, pch = c(21, 8, 21, 8, 21, 16, 21, 16, 21, 16, 21, 16, 21, 16, 21, 16, 21, 16), yaxt = "n", ylab = "", main = "Women", xlim = c(-.50, .50), ylim = c(.5, 35.3), cex.axis = .8, cex.main = 1)
abline(v = 0, lty = 2)








#####REGRESSIONS FOR CANDIDATE VOTE CHOICE BY RACE (OA-L)#####

#Organize Data

smallstackedb <- subset(stacked.b, select = c(id, candvote, treatment, candgender, candrace, candideology, candpolicy, female), treatment == 0) #Create a subset that includes only the variables needed for the regressions
head(smallstackedb)
nrow(smallstackedb)

smallstackedb$candrace <- as.factor(smallstackedb$candrace) #Create factor variables, white, moderate, good jobs, male as the baseline
smallstackedb$candpolicy <- as.factor(smallstackedb$candpolicy)
smallstackedb$candideology <- as.factor(smallstackedb$candideology)

women.black <- subset(smallstackedb, female == 1)
men.black <- subset(smallstackedb, female == 0)

smallstackedw <- subset(stacked.w, select = c(id, candvote, treatment, candgender, candrace, candideology, candpolicy, female), treatment == 0) #Create a subset that includes only the variables needed for the regressions
head(smallstackedw)
nrow(smallstackedw)

smallstackedw$candrace <- as.factor(smallstackedw$candrace) #Create factor variables, white, moderate, good jobs, male as the baseline
smallstackedw$candpolicy <- as.factor(smallstackedw$candpolicy)
smallstackedw$candideology <- as.factor(smallstackedw$candideology)

women.white <- subset(smallstackedw, female == 1)
men.white <- subset(smallstackedw, female == 0)

women.black.control <- subset(women.black, treatment == 0)
women.black.treat <- subset(women.black, treatment == 1)
men.black.control <- subset(men.black, treatment == 0)
men.black.treat <- subset(men.black, treatment == 1)

women.white.control <- subset(women.white, treatment == 0)
women.white.treat <- subset(women.white, treatment == 1)
men.white.control <- subset(men.white, treatment == 0)
men.white.treat <- subset(men.white, treatment == 1)

mwc <- zelig(candvote ~ candgender + candrace + candideology + candpolicy, model = "logit.gee", id = "id", data = men.white.control, cite = FALSE)
summary(mwc)
nrow(na.omit(men.white.control))

wwc <- zelig(candvote ~ candgender + candrace + candideology + candpolicy, model = "logit.gee", id = "id", data = women.white.control, cite = FALSE)
summary(wwc)
nrow(na.omit(women.white.control))

mbc <- zelig(candvote ~ candgender + candrace + candideology + candpolicy, model = "logit.gee", id = "id", data = men.black.control, cite = FALSE)
summary(mbc)
nrow(na.omit(men.black.control))

wbc <- zelig(candvote ~ candgender + candrace + candideology + candpolicy, model = "logit.gee", id = "id", data = women.black.control, cite = FALSE)
summary(wbc)
nrow(na.omit(women.black.control))



#####REGRESSIONS FOR CANDIDATE VOTE CHOICE BY RACE AND GENDER BASED ON TREATMENT ACCEPTANCE (OA-M)#####

#Organize Data

treated.stack.white <- subset(stacked.w, treatment == 1)
treated.stack.white

treated.stack.black <- subset(stacked.b, treatment == 1)
treated.stack.black

smallstackedb <- subset(treated.stack.black, select = c(id, candvote, clintonid, candgender, candrace, candideology, candpolicy, female)) #Create a subset that includes only the variables needed for the regressions
head(smallstackedb)
nrow(smallstackedb)

smallstackedb$candrace <- as.factor(smallstackedb$candrace) #Create factor variables, white, moderate, good jobs, male as the baseline
smallstackedb$candpolicy <- as.factor(smallstackedb$candpolicy)
smallstackedb$candideology <- as.factor(smallstackedb$candideology)

women.black <- na.omit(subset(smallstackedb, female == 1))
men.black <- subset(smallstackedb, female == 0)

smallstackedw <- subset(treated.stack.white, select = c(id, candvote, clintonid, candgender, candrace, candideology, candpolicy, female)) #Create a subset that includes only the variables needed for the regressions
head(smallstackedw)
nrow(smallstackedw)

smallstackedw$candrace <- as.factor(smallstackedw$candrace) #Create factor variables, white, moderate, good jobs, male as the baseline
smallstackedw$candpolicy <- as.factor(smallstackedw$candpolicy)
smallstackedw$candideology <- as.factor(smallstackedw$candideology)

women.white <- subset(smallstackedw, female == 1)
men.white <- subset(smallstackedw, female == 0)



m1b <- zelig(candvote ~ candgender*clintonid + candrace*clintonid + candideology*clintonid + candpolicy*clintonid, model = "logit.gee", id = "id", data = men.white, cite = FALSE)
summary(m1b)
nrow(na.omit(men.white))
count(unique(men.white$id))

m2b <- zelig(candvote ~ candgender*clintonid + candrace*clintonid + candideology*clintonid + candpolicy*clintonid, model = "logit.gee", id = "id", data = women.white, cite = FALSE)
summary(m2b)
nrow(na.omit(women.white))
count(unique(women.white$id))

m3b <- zelig(candvote ~ candgender*clintonid + candrace*clintonid + candideology*clintonid + candpolicy*clintonid, model = "logit.gee", id = "id", data = men.black, cite = FALSE)
summary(m3b)
nrow(na.omit(men.black))
count(unique(men.black$id))

m3b <- zelig(candvote ~ candgender*clintonid + candrace*clintonid + candideology*clintonid + candpolicy*clintonid, model = "logit.gee", id = "id", data = women.black, cite = FALSE)
summary(m3b)
nrow(na.omit(women.black))
count(unique(women.black$id))


#####REGRESSIONS FOR CANDIDATE VOTE CHOICE BY RACE AND GENDER BASED ON JUST ACCEPTANCE (OA-N)#####

#Organize Data

smallstackedb <- subset(stacked.b, select = c(id, candvote, clintonid, candgender, candrace, candideology, candpolicy, female)) #Create a subset that includes only the variables needed for the regressions
head(smallstackedb)
nrow(smallstackedb)

smallstackedb$candrace <- as.factor(smallstackedb$candrace) #Create factor variables, white, moderate, good jobs, male as the baseline
smallstackedb$candpolicy <- as.factor(smallstackedb$candpolicy)
smallstackedb$candideology <- as.factor(smallstackedb$candideology)

women.black <- subset(smallstackedb, female == 1)
men.black <- subset(smallstackedb, female == 0)

smallstackedw <- subset(stacked.w, select = c(id, candvote, clintonid, candgender, candrace, candideology, candpolicy, female)) #Create a subset that includes only the variables needed for the regressions
head(smallstackedw)
nrow(smallstackedw)

smallstackedw$candrace <- as.factor(smallstackedw$candrace) #Create factor variables, white, moderate, good jobs, male as the baseline
smallstackedw$candpolicy <- as.factor(smallstackedw$candpolicy)
smallstackedw$candideology <- as.factor(smallstackedw$candideology)

women.white <- subset(smallstackedw, female == 1)
men.white <- subset(smallstackedw, female == 0)





m1b <- zelig(candvote ~ candgender*clintonid + candrace*clintonid + candideology*clintonid + candpolicy*clintonid, model = "logit.gee", id = "id", data = men.white, cite = FALSE)
summary(m1b)
nrow(na.omit(men.white))

m2b <- zelig(candvote ~ candgender*clintonid + candrace*clintonid + candideology*clintonid + candpolicy*clintonid, model = "logit.gee", id = "id", data = women.white, cite = FALSE)
summary(m2b)
nrow(na.omit(women.white))

m3b <- zelig(candvote ~ candgender*clintonid + candrace*clintonid + candideology*clintonid + candpolicy*clintonid, model = "logit.gee", id = "id", data = men.black, cite = FALSE)
summary(m3b)
nrow(na.omit(men.black))

m4b <- zelig(candvote ~ candgender*clintonid + candrace*clintonid + candideology*clintonid + candpolicy*clintonid, model = "logit.gee", id = "id", data = women.black, cite = FALSE)
summary(m4b)
nrow(na.omit(men.black))










